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ABSTRACT 


An  imaging  radar,  like  ISAR,  offers  a  combatant  the  capability  to  perform 
long  range  surveillance  with  high  quality  imagery  for  positive  target 
identification.  Extending  this  attractive  feature  to  the  battle  damage  assessment 
problem  (BDA)  gives  the  operator  instant  viewing  of  the  target’s  behavior  when 
it  is  hit.  As  a  consequence,  immediate  and  decisive  action  can  be  quickly  taken 
(if  required).  However,  the  conventional  Fourier  processing  adopted  by  most 
ISAR  systems  does  not  provide  adequate  time  resolution  to  capture  the  target’s 
dynamic  responses  during  the  hit.  As  a  result,  the  radar  image  becomes 
distorted.  To  improve  the  time  resolution,  time-frequency  transform  (TFT) 
methods  of  ISAR  imaging  have  been  proposed.  Unlike  traditional  Fourier- 
based  processing,  TFT’s  allows  variable  time  resolution  of  the  entire  event  that 
falls  within  the  ISAR  coherent  integration  period  to  be  extracted  as  part  of  the 
imaging  process.  We  have  shown  in  this  thesis  that  the  use  of  linear  Short 
Time-Frequency  Transforms  allows  the  translational  response  of  the  aircraft 
caused  by  a  blast  force  to  be  clearly  extracted.  The  TFT  extracted  images  not 
only  tell  us  how  the  aircraft  responds  to  a  blast  effect  but  also  provides 
additional  information  about  the  cause  of  image  distortion  in  the  traditional 
ISAR  display. 


v 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


VI 


TABLE  OF  CONTENTS 


I.  INTRODUCTION . 1 

A.  ISAR  FOR  BATTLE  DAMAGE  ASSESSMENT  (BDA) . 1 

B.  RADAR  TARGET  IMAGING  USING  ISAR . 1 

C.  RESPONSE  OF  TARGET  BEING  HIT . 3 

D.  JOINT  TIME-FREQUECY  METHOD  OF  ISAR  IMAGING  FOR 

BDA . 3 

E.  THESIS  OUTLINE . 5 

II.  INVERSE  SYNTHETIC  APERTURE  RADAR  (ISAR)  THEORY  AND 

SIGNAL  PROCESSING . 7 

A.  RADAR  IMAGING  AND  ITS  ASSOCIATED  RESOLUTIONS . 7 

B.  ISAR  SYSTEM  ARCHITECTURE . 8 

C.  THEORY  OF  MOVING  TARGET  IMAGING  IN  ISAR . 10 

1.  ISAR  Geometry  and  its  Relationship  with  the  Returned 

Signal . 10 

2.  Signal  Representation . 12 

3.  ISAR  Data  Collection  and  Processing . 13 

D.  ACHIEVING  HIGH  RANGE  RESOLUTION  THROUGH 

STEPPED  FREQUENCY  WAVEFORM  (SFW)  SIGNAL 
PROCESSING . 15 

1.  Synthetic  Range  Profile  Processing . 15 

2.  Unambiguous  Range  and  Range  Resolution . 18 

3.  Effect  of  Target  Velocity . 18 

4.  Range  Offset  and  Range  Walk . 19 

E.  ACHIEVING  HIGH  CROSS-RANGE  RESOLUTION  USING 

SYNTHETIC  APERTURE . 21 

1 .  ISAR  Theory  from  an  Aperture  Viewpoint . 21 

2.  Range  Doppler  Imaging  and  Cross  Resolution . 23 

3.  Doppler  Frequency  Shift  Due  to  Target  Motion . 25 

4.  Distortion  Produced  by  Target  Rotation . 26 

III.  BLAST  LOADING  ON  AIRCRAFT . 29 

A.  BLAST  IMPACT  ON  AIRBORNE  TARGETS . 29 

B.  CHARACTERISTIC  OF  EXPLOSIVE  BLAST  IN  AIR . 30 

1 .  Formation  of  the  Blast  Wave . 30 

2.  Characteristic  of  the  Blast  Wave . 32 

3.  Peak  Overpressure . 34 

4.  Arrival  Time . 34 

5.  Blast  Duration . 35 

6.  Impulse  Per  Unit  Area . 36 

7.  Hopkinson  Scaling  Law . 37 

vii 


C.  DYNAMIC  LOADING  OF  BLAST  WAVE  ON  AIRCRAFT . 38 

1.  Translational  and  Rotational  Responses . 38 

2.  Approach  to  Numerical  Simulation  for  ISAR  Imaging . 42 

D.  BLAST  WAVE  EFFECT  ON  AIRCRAFT  AND  ISAR  IMAGING . 44 

IV.  TIME-FREQUENCY  TRANSFORM  METHODS  OF  ISAR  IMAGING . 45 

A.  TIME-FREQUENCY  REPRESENTATIONS  OF  SIGNALS . 45 

B.  SHORT  TIME-FREQUENCY  TRANSFORM  (STFT) . 46 

1.  Definition . 46 

2.  Properties  of  STFT . 46 

3.  Time  and  Frequency  Resolution . 47 

4.  Shape  of  the  Window  Function . 50 

C.  BI-LINEAR  TIME-FREQUENCY  TRANSFORM . 50 

D.  TIME-FREQUENCY  BASED  IMAGE  FORMATION . 52 

V.  DEMONSTRATION  OF  BLAST  EFFECT  CAPTURINGBY  ISAR  USING 

TIME-FREQUENCY  TRANSFORM . 57 

A.  BLAST  EFFECTS  AND  THE  ISAR  SIMULATION  MODEL . 57 

B.  SIMULATION  OF  AIRCRAFT  RESPONSE  TO  BLAST  WAVE 

INTERACTION . 58 

1.  Blast  Pressure  Distribution  on  Aircraft . 58 

2.  Geometry  of  the  Problem  and  its  Analytical  Solution . 60 

3.  Aircraft  Profile . 62 

4.  Blast  Parameter . 63 

C.  SIMULATION  OF  ISAR  IMAGING  USING  TIME-FREQUENCY 

TRANSFORM . 64 

1.  Target  Model  Representation . 64 

2.  Radar  Signal  Generation . 64 

3.  ISAR  Image  Processor . 66 

VI.  RESULTS  AND  DISCUSSIONS . 67 

A.  TIME-FREQUENCY  REPRESENTATION  OF  ISAR  IMAGE . 67 

B.  TRANSLATIONAL  RESPONSE  OF  AIRCRAFT  DUE  TO  BLAST 

EFFECT . 69 

1 .  Blast  Characteristic . 69 

2.  ISAR  Image  from  Conventional  FFT . 70 

3.  ISAR  image  from  STFT . 71 

C.  SELECTION  OF  WINDOW  FUNCTIONS . 76 

1.  Resolution  of  Rectangular  and  Gaussian  Window 

Functions . 76 

2.  Effect  of  Using  Rectangular  Window  Function  in  STFT 

Imaging . 77 

3.  Effect  of  Using  Gaussian  Window  Function  in  STFT 

Imaging . 79 

D.  EFFECT  OF  RADIAL  VELOCITY  VARIATION  ON  ISAR  IMAGE . 82 


viii 


VII.  CONCLUSION  AND  RECOMMENDATIONS . 87 

A.  CONCLUSION . 87 

1.  Approach  to  Model  the  Aircraft  Response  When  Hit . 88 

2.  Approach  to  Overcome  the  Image  Distortion . 89 

B.  RECOMMENDATIONS  FOR  FURTHER  WORK . 90 

1 .  Complete  Model  of  the  Blast  Loading  on  the  Aircraft . 90 

2.  Study  of  Target  Break  Up  Process . 91 

3.  Improving  the  Time  and  Frequency  Resolution  of  TFT 91 

APPENDIX  A:  MATLAB  CODES . 93 

LIST  OF  REFERENCES . 97 

INITIAL  DISTRIBUTION  LIST . 99 


IX 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


x 


LIST  OF  FIGURES 


Figure  1  A  simple  illustration  of  SWF  ISAR  architecture . 9 

Figure  2  ISAR  geometry  in  two-dimensional  space . 10 

Figure  3  Row-column  decomposition  for  synthetic  ISAR  imaging . 14 

Figure  4  SFW  Pulse  and  echo  pulses . 16 

Figure  5  SAR  and  ISAR  equivalence . 22 

Figure  6  Radial  velocity  from  a  point  scatterer  on  a  rotating  target . 24 

Figure  7  Development  of  explosive  shock,  (a)  Initial  pressure  pulse,  (b) 

successive  configurations,  and  (c)  Final  form  when  fully 

developed . 31 

Figure  8  Typical  pressure-distance  curves  for  successive  times  after  an 

explosion . 32 

Figure  9  Typical  pressure-time  curves  for  explosive  blast  wave . 32 

Figure  10  3D  geometry  of  the  interaction  of  blast  wave  on  the  aircraft 

surface . 39 

Figure  1 1  3D  geometry  of  the  ISAR  staring  at  aircraft . 42 

Figure  12  Approach  to  generate  of  target  motion  for  ISAR  simulation . 43 

Figure  13  Example  of  a  time-varying  signal  with  changing  frequency . 48 

Figure  14  Effect  of  Time  and  frequency  resolution . 49 

Figure  15  ISAR  image  processing  based  of  TFT  method . 54 

Figure  16  Extraction  of  ISAR  Image  from  the  Doppler-time-Range  Matrix . 55 

Figure  17  2D  Geometry  between  ISAR,  aircraft  and  blast  point . 60 

Figure  18  Aircraft  area  profile . 63 

Figure  19  Target  model  representation  and  its  ISAR  Image  (from  FFT 

process) . 67 

Figure  20  Sample  of  ISAR  from  STFT . 68 

Figure  21  Target  acceleration  and  velocity  variation  (Blown  up  version) . 69 

Figure  22  Target  range  variation . 69 

Figure  23  ISAR  image . 70 


XI 


Figure  24  STFT  ISAR  images  for  observing  the  blast  effect . 72 

Figure  25  Acceleration  and  Velocity  variation  for  2nd  case . 73 

Figure  26  ISAR  image  from  FFT  for  2nd  Case . 74 

Figure  27  STFT  ISAR  images  for  2nd  case . 75 

Figure  28  STFT  ISAR  Time-history  images  with  width  set  to  2  m-cell . 78 

Figure  29  STFT  ISAR  Time-history  images  with  width  set  to  55  m-cell . 79 

Figure  30  Time  and  frequency  resolution  comparison  chart  for  rectangular 

and  Gaussian  window  functions . 80 

Figure  31  STFT  ISAR  Time-history  images  with  Gaussian  window  function . 81 

Figure  32  Extent  of  cell  migration  due  to  velocity . 82 


xii 


LIST  OF  TABLES 


Table  1 .  Radar  Parameters 


65 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


XIV 


ACKNOWLEDGMENTS 


The  author  would  like  to  acknowledge  Professor  Brett  Borden  for  his 
guidance,  patience  and  motivation  throughout  the  thesis  process  to  keep  the 
work  on  the  right  track.  In  addition,  the  author  would  like  to  thank  Professor 
Ronald  Brown  for  sharing  his  valuable  expert  knowledge  in  the  area  of 
explosive  study.  Also,  the  author  would  like  to  give  credit  to  Professor  Donald 
Walters  for  providing  critical  comments  to  help  in  improving  the  thesis  work  and 
in  the  report.  Without  them,  greater  understanding  of  the  subjects  involved  in 
this  thesis  work  could  not  be  achieved. 


xv 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


XVI 


I.  INTRODUCTION 


A.  ISAR  FOR  BATTLE  DAMAGE  ASSESSMENT  (BDA) 

A  radar  target1  image  is  generated  from  the  reflectivity  data  collected  by  a 
radar  platform  as  the  target  is  observed  from  a  set  of  viewing  angles.  For 
Synthetic  Aperture  Radar  (or  SAR),  the  radar  platform  moves  to  give  the  different 
target  viewing  angles,  whereas  Inverse  Synthetic  Aperture  Radar  (ISAR)  makes 
use  of  the  target’s  angular  rotations  such  as  roll,  pitch  and  yaw  to  collect  the 
viewing  angles  for  image  generation. 

ISAR  imaging  opens  up  the  new  possibility  of  using  radar  measurements 
to  construct  high  quality  images  for  target  recognition  purposes  while  the  target 
(usually  aircraft  and  ships)  is  moving  and  the  radar  platform  remains  essentially 
stationary2.  Operationally,  ISAR  offers  enhanced  target  recognition  at  ranges 
that  are  often  not  achievable  by  electro-optic  systems  or  Forward  Looking 
Infrared  (FUR)  alone.  An  extension  of  this  long-range  target  recognition 
capability  of  the  ISAR,  apart  from  deploying  for  surveillance  purpose,  is  to  use  it 
to  perform  Battle  Damage  Assessment  (BDA). 

The  significant  advantage  of  using  ISAR  for  BDA  (over  other  methods)  is 
its  ability  to  provide  the  operator  instantaneous  viewing,  in  highly  resolved  detail, 
of  the  behavior  of  the  target  as  a  kinetic  energy  weapon  hits  it.  This  gives  the 
operator  immediate  evaluation  of  the  extent  of  the  damage  caused  by  the  hit. 
Follow-up  counter  actions  can  be  quickly  taken  if  needed  (all  being  performed 
beyond  the  target’s  weapon  range). 


B.  RADAR  TARGET  IMAGING  USING  ISAR 

To  construct  a  radar  image  from  the  reflectivity  data,  most  ISAR  systems 
process  the  data  by  performing  a  two  dimensional  (2D)  Fast  Fourier  Transform 

1  The  target  can  be  any  object  of  interest  to  the  radar  operator. 

2  The  platform  can  be  in  motion  but  its  motion  has  to  be  accurately  measured  (usually  by  the 
combination  of  Inertial  Navigation  System  (INS)  and  Global  Positioning  System  (GPS),  and 
compensated  accordingly. 
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(FFT);  the  first  FFT  generates  the  synthetic  range  profile  of  the  target,  and  the 
second  FFT  extracts  the  Doppler  information  for  each  range  cell.  Ideally,  the 
resultant  2D  Range-Doppler  plot  following  this  process  would  closely  resemble 
an  outline  of  the  target,  arising  from  the  various  scattering  points  around  the 
target.  The  quality  of  the  generated  ISAR  image  is  directly  related  to  the 
resolution  imposed  by  the  processing  system.  High  range  resolution  can  be 
achieved  through  pulse  compression,  and  Stepped  Frequency  Waveforms 
(SFW)  are  the  most  common  methods  used  today.  High  cross-range  resolution 
depends  strongly  on  the  variation  in  viewing  angles  (equivalently  Doppler 
resolution)  which,  in  turn,  increases  with  higher  angular  rate  of  target  rotation  and 
target  dwell  time. 

However,  a  subtlety  that  affects  the  quality  of  the  ISAR  image  results  from 
the  requirement  of  minimum  variation  in  the  range  over  the  different  viewing 
angles  between  the  radar  platform  and  the  target  for  cross-range  synthesis. 
Significant  variation  in  range  (due  to  the  presence  of  translational  velocity 
components,  non-linear  or  large  variations  in  angular  velocity,  or  combination  of 
both)  will  introduce  an  additional  quadratic  phase  error  that  tends  to  distort  or  blur 
(known  also  as  the  defocusing  effect)  the  image  during  the  signal  processing 
stage.  This  phenomenon  is  called  “Range  Walk”  by  the  SAR  imaging 
community,  and  it  acts  to  impose  an  upper  limit  on  the  effective  observed  target 
rotation  interval  and — because  this  interval  is  determined  by  rotation  rate — range 
walk  also  bounds  the  integration  time. 

Normally,  if  the  radial  velocity  remains  essentially  constant,  the  motion  of 
the  target  can  be  accurately  estimated  (and  compensated  for)  to  reduce  the 
effects  of  range  walk.  The  angular  velocity,  on  the  other  hand,  is  needed  for 
ISAR  imaging  and  is  generally  assumed  to  be  constant  and  reasonably  small 
over  the  entire  data  integration  time.  In  most  (aircraft  or  ship)  target’s  roll,  pitch 
and  yaw  behaviors,  this  is  approximately  valid  for  the  observation  period.  The 
converse  is  not  true  though  if  the  target  experiences  any  acceleration. 
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C.  RESPONSE  OF  TARGET  BEING  HIT 

While  the  target  is  being  hit,  the  blast  force  from  the  explosion  of  the 
weapon — or  weapon  impact  itself — acting  on  the  target  structure  will  cause  it  to 
undergo  translational  and  rotational  motions,  as  well  as  structural  vibration.  The 
time  duration  of  this  force  acting  on  the  target  is  often  very  short  (on  order  of 
milliseconds)  but  the  impulse  is  large  enough  to  cause  significant  alterations  in 
both  radial  and  angular  velocity  of  the  target.  As  a  consequence  of  this  short 
impulse  period,  the  radar  may  not  have  adequate  time  to  compensate  for  these 
changes  in  velocities  and  hence  errors  are  introduced  into  the  image  processing. 
These  errors  could  account  for  the  distortion  in  most  ISAR  image  observed  when 
the  target  is  being  hit. 

In  addition  to  the  translational,  rotation  and  vibration  effects,  the  sharp  and 
high  velocity  fragments  from  the  weapon  as  it  explodes,  or  possible  a  Kinetic 
Energy  (KE)  round  itself,  are  likely  to  penetrate  the  target  causing  it  to  break  up. 
The  smaller  broken  target  pieces  now  begin  to  scatter  and  reflect 
electromagnetic  signals  back  to  the  radar  as  they  move  at  much  higher  velocities 
(from  the  energy  acquired)  than  the  target  (which  has  larger  mass).  The  returned 
scattered  EM  signal  from  these  pieces  can  behave  like  point  scatterers,  dihedral, 
trihedral  or  multi-facet  corner  reflectors,  or  specular  reflectors  depending  on  the 
dynamic  nature  of  the  break-up  process.  It  is  difficult  to  predict  target  structure 
changes  without  an  in-depth  analysis  based  on  explosive  and  structural  impact 
studies.  In  most  cases,  numerical  simulations  such  as  finite  element  analysis  are 
often  used  in  these  studies  and  the  results  obtained  are  only  specific  to  the 
particular  scenario  defined.  Due  to  time  constraints,  this  topic  will  not  be  pursued 
in  the  thesis  work. 

D.  JOINT  TIME-FREQUECY  METHOD  OF  ISAR  IMAGING  FOR  BDA 

Time-frequency  methods  for  signal  analysis  (also  known  as  the  Time- 
Frequency  Transform  or  TFT)  are  similar  to  conventional  Fourier  Transform  (FT) 
schemes  except  that  the  FT  is  taken  after  the  auto-correlation  of  the  signal  with  a 
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weighted  window  function.  The  advantage  of  this  analysis  method  is  that  it  allows 
decomposition  of  the  frequency  spectrum  of  time-varying  signals  in  shorter 
window  time  frames,  and  provides  an  added  dimension  for  examining  the 
dynamic  behavior  of  the  signal  as  it  varies  over  time.  The  size  and  shape  of  the 
weighting  window  function  can  be  altered  to  fit  the  specific  analysis  requirement 
for  the  signal. 

Time-frequency  methods  in  signal  analysis  have  aroused  much  interest  in 
the  signal  processing  community.  The  initial  application  was  in  spectrogram 
analysis  of  complex  human  speech.  However,  in  recent  years,  further  extension 
of  this  time-frequency  method  to  ISAR  was  actively  advocated  by  Chen  [2]  and 
Jae  and  Thomas  and  Flores  [3]  —  particularly  in  the  area  of  image  generation, 
motion  compensation  and  micro-Doppler  target  vibration  studies.  In  image 
generation,  Chen  [2]  has  reportedly  used  the  TFT  to  construct  ISAR  image  of 
targets  with  high  rotation  rate  by  helping  to  overcome  the  migration  of  individual 
scattering  points  from  one  range  cell  to  another.  Both  Jae,  Thomas  and  Flores 
[3]  and  Chen  [2]  have  also  proposed  the  use  of  TFTs  to  resolve  multiple  targets 
appearing  in  ISAR  images. 

So  far,  the  possibility  of  using  the  TFT  to  examine  the  dynamic  behavior  of 
a  target  in  an  ISAR  image  as  the  target  is  being  acted  upon  by  an  impulse  blast 
force  has  yet  to  be  explored  by  other  TFT-centered  research  efforts.  If  the  target 
is  assumed  to  be  a  rigid  body,  it  would  experience  sudden  translational  and 
rotational  motion  variation  under  blast  wave  loading.  These  motions  will  be 
measured  during  the  data  collection  process  if  the  radar  is  continuously  staring  at 
the  target.  This  thesis  aims  to  apply  the  appropriate  TFT  window  function  to  the 
collected  ISAR  data  and  demonstrate  the  possibility  of  extracting  these  target 
motions  using  TFT  methods  of  image  construction. 


4 


E.  THESIS  OUTLINE 

The  thesis  is  organized  in  the  following  manner: 

Chapter  II  serves  to  develop  the  theory  of  ISAR  image  processing.  This 
chapter  will  touch  on  how  high  quality  image  resolution  is  achieved,  and  the 
conditions  necessary  to  maintain  the  quality  of  the  image  generated. 

In  Chapter  III,  the  effects  of  blast  waves  interacting  with  a  target  are 
described  with  the  goal  of  modeling  and  simulation.  The  subject  target  for  this 
study  is  an  aircraft  which  is  treated  as  a  rigid  body  to  simplify  the  analysis.  The 
translational  and  rotational  motions  resulting  from  force  interaction  are  examined. 

Chapter  IV  introduces  the  Time-Frequency  Transform  (TFT)  Method  of 
analysis,  and  how  it  can  be  applied  to  ISAR  signal  processing  and  image 
generation.  The  discussion  in  this  chapter  is  largely  based  on  the  work  of  Chen 
[2]  in  TFT  method  for  ISAR  imaging  and  Jae  and  Thomas  and  Flores  [3]  for 
motion  compensation. 

To  demonstrate  the  feasibility  of  using  TFT  methods  for  ISAR  imaging, 
and  their  advantage  over  conventional  FFT  based  methods  for  BDA  (oweing  to 
the  dynamic  behavior  of  the  aircraft  as  it  is  hit),  a  simple  simulation  program  is 
developed.  Chapter  V  describes  how  this  program  is  put  together  for  modeling 
the  radar  target  returns  from  an  aircraft,  the  effect  of  the  blast  loading  on  the 
aircraft  and  the  construction  of  ISAR  images  using  TFT  concepts.  This  chapter 
also  discusses  the  parameters  and  conditions  used  in  the  simulation. 

Chapter  VI  follows  on  to  discuss  the  results  of  the  simulation  and  its 
significance. 

Chapter  VII  is  the  concluding  chapter  and  summarizes  the  major  findings 
of  the  thesis  work  and  provides  recommendations  for  similar  future  work. 
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II.  INVERSE  SYNTHETIC  APERTURE  RADAR  (ISAR) 
THEORY  AND  SIGNAL  PROCESSING 

A.  RADAR  IMAGING  AND  ITS  ASSOCIATED  RESOLUTIONS 

Radar  has  long  been  recognized  for  its  ability  to  “detect,  locate  and 
identify  targets  at  great  distance  and  in  all  kind  of  weather  conditions”  [1]  and,  to 
date,  it  is  hardly  matched  by  other  sensors.  Compared  to  Infrared  (IR),  visible  or 
Ultra-violet  sensors,  the  longer  wavelengths  (by  four  to  six  orders  of  magnitude) 
allow  the  radar  energy  to  propagate  through  most  of  the  atmosphere  with  little 
loss,  even  in  poor  weather  conditions.  However,  the  long  wavelengths  limit  the 
radar’s  imaging  resolution,  which  is  proportional  to  the  ratio  of  the  wavelength  to 
the  aperture  size,  and  hence  affects  the  image  quality.  To  achieve  the  fine 
image  resolution  in  optics,  a  lens  system  with  aperture  size  of  up  to  several 
centimeters  is  adequate.  But  a  radar  system  would  require  an  antenna  aperture 
diameter  to  stretch  tens  of  kilometers  to  attain  the  same  resolution.  Thus, 
alternate  methods  to  accomplish  this  must  be  found. 

Inverse  Synthetic  Aperture  Radar  (ISAR)  combines  the  use  of  pulse 
compression,  high  Pulse  Repetition  Frequency  (PRF)  and  the  target  motions 
(particularly  its  angular  rotation)  to  generate  high  range  and  cross-range 
resolution  target  images.  Range  resolution  is  associated  with  the  width  of  the 
target  pulse  in  the  time  domain.  Hence,  narrower  pulses  will  have  better 
resolution.  But  such  radar  systems  require  higher  peak  power  to  ensure  enough 
energy  returned  from  the  target  for  signal  processing.  To  mitigate  this,  pulse 
compression  is  mostly  employed  and  for  ISAR  a  common  technique  uses 
Stepped  Frequency  Waveforms  (SFW).  High  cross-range  resolution  can  be 
achieved  by  observing  the  phase  variations  of  the  signals  returned  from  different 
target  viewing  angles  over  a  coherent  processing  or  integration  period.  It  can  be 
shown,  under  rather  stringent  conditions  (in  subsequent  sections),  that  the 
changes  in  phase  due  the  rotation  of  the  target  can  be  mapped  in  one-to-one 
correspondence  to  the  cross-range  position  of  the  target  scattering  points. 
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B.  ISAR  SYSTEM  ARCHITECTURE 

Figure  1  below  shows  the  architecture  of  a  simple  SFW  ISAR  system.  At 
the  RF  stage,  a  series  of  high  PRF  narrow  pulses  is  used  to  modulate  the  RF 
carrier  that  varies  in  frequency  steps  of  Af  at  every  pulse  interval,  T2,  which  is 
reciprocal  to  the  PRF.  A  burst  of  transmitted  pulses  consists  of  N  stepped 
frequency  pulses  varying  from  £  to  fc  +  (n-l)Af,  where  n  =  1,  2,  N..  For  ISAR 

image  processing,  M  bursts  are  required  to  form  an  image.  The  size  of  the 
image  is  thus  determined  by  the  M  x  N  sequence  of  the  pulses,  and  the 
integration  time,  T,  needed  to  construct  the  image  (if  we  ignored  the  processing 
delays)  is  then  T  =  M  x  N  x  T2.  The  MxN  SFW  pulse  sequence  represents  the 
Doppler-Range  mapping  for  the  image. 


f  -  Intermediate  frequency 
fc  -  Carrier  Frequency 
Af  -  Stepped  frequency 
n  =  1,  2,  ....  N 


(a)  RF  transmission  and  detection  stages  (After:  Figure  5.2(b),  pp  202,  Wehner  [4]) 
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Complex  I  and  Q  samples 


M  Time  Histories 


(b)  Signal  and  image  processing  stages  (After:  Figure  5.2,  pp  100,  Chen  [2]) 


Figure  1  A  simple  illustration  of  SWF  ISAR  architecture. 

The  returned  echo  signal  captured  by  the  receiver  immediately  after  each 
pulse  transmission  is  down-converted  to  an  intermediate  frequency  (IF)  and 
passed  along  for  quadrature  detection.  The  quadrature  detector  detects  the 
baseband  signal  in  the  form  of  in-phase  (I)  and  quadrature  phase  (Q)  signals. 
These  signals  are  then  sampled  and  held  as  an  M  x  N  data  array  for  signal 
processing.  The  sampled  output  data  represents  the  time-history  of  the  target’s 
reflectivity  in  the  frequency  domain. 

At  the  signal  processing  stage  shown  in  Figure  1(b),  corrections  for  radar 
system  “phase  and  amplitude  ripples”  (Wehner  [4])  are  first  applied,  followed  by 
range  alignment,  translational  and  rotational  motion  compensations  in  the  form  of 
range  and  Doppler  tracking.  Motion  compensations  are  employed  to  minimize 
the  induced  quadratic  phase  error  and  hence  “Range  Walk”  effect.  The  actual 
image  construction  begins  with  the  synthetic  range  profile  processing  which  is 
simply  formed  by  taking  an  inverse  Fourier  transform  of  the  N  pulses  for  each  m 
burst  where  m  =  1,  2,  ....  M.  Cross-range  information  is  then  extracted  by  another 
Fourier  process  -  a  direct  Fourier  transform  this  time  -  along  the  column  of  M 
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bursts  for  each  n  range  cell  index.  In  practice,  the  synthetic  range  profile  and 
cross-range  (Doppler  history)  information  are  simultaneously  processed  by  taking 
a  two-dimensional  Fast  Fourier  Transform3  (FFT)  of  the  compensated  MxN data 
array  to  lessen  the  computation  load. 


C.  THEORY  OF  MOVING  TARGET  IMAGING  IN  ISAR 

1.  ISAR  Geometry  and  its  Relationship  with  the  Returned  Signal 

For  illustration  purposes,  consider  a  simplified  two-dimensional  (2D) 
geometry  between  the  radar  and  an  airborne  target  as  shown  below  in  Figure  2. 
The  aircraft  is  assumed  to  have  translational  and  rotational  motion  only  in  the 
two-dimensional  plane  relative  to  a  stationary  radar  platform.  It  can  be  easily 
extended  to  the  actual  three-dimensional  case  to  include  the  effect  of  the  motions 
along  and  about  the  other  axes  with  some  careful  considerations. 


co  rotation 


Figure  2  ISAR  geometry  in  two-dimensional  space. 

(After:  Figure  5.1,  pp95,  Chen  [2]) 

3  It  can  be  shown  that  the  process  of  IDFT  and  subsequently  DFT  is  equivalent  to  taking  a 
two-dimensional  FFT  with  the  exception  of  a  phase  factor  in  front.  In  most  radar  signal 
processing,  we  can  ignore  this  phase  factor  since  our  interest  is  largely  on  the  power  spectral 
density  or  absolute  intensity  of  the  returns. 
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Following  Chen  [2],  we  designate  the  radar  to  be  located  at  the  origin  of 
the  U-V  coordinate  system,  known  as  radar  coordinates.  The  target  can  be 
described  by  a  collection  of  scattering  points  (or,  “scattering  centers”)  in  its  own 
reference  frame,  x-y,  known  as  the  target  coordinates.  The  center  of  rotation  of 
the  target  is  located  at  the  origin  of  this  x-y  coordinate  system.  To  describe  the 
rotation  of  the  target  independent  of  the  target  scattering  center  distribution  in  x-y 
coordinates,  another  reference  coordinate  system  denoted  as  X-Y  is  introduced 
with  its  origin  coincident  with  the  target  coordinates.  The  two  axes  in  X-Y  are 
parallel  to  U  and  V  axes  respectively  so  that  a  pure  translational  transformation 
can  be  used  to  describe  their  relationship.  Then,  it  follows  that  the 
transformation  of  x-y  reference  to  X-Y  coordinate  can  be  described  by  a  pure 
rotation.  Also,  to  simplify  the  illustration  further,  we  assume  a  =  0  (the  radar 
target  line  of  sight  (LOS)  is  along  the  U-axis  of  the  radar  reference). 

Assume  that  at  t  =  0  the  target  range  (the  distance  from  the  radar  to  the 
center  of  rotation  of  the  target)  is  R.  Hence,  the  distance  between  the  radar  and  a 
scattering  point,  P(x,y),  can  be  expressed  as 


RP  = 


[R  +  x cos  90  -  vsin  90 )  +  ( y  cos  90  +  xsin  9(j ) 


1/2 


=  [i?2  +x 2  +  y°  +  2R(xcos90  -  jsin^)] 
=  R  +  x  cos  90  -  v  sin  90 


(2.1) 


where  9a  is  the  initial  target  rotation  angle. 

If  the  target  possesses  both  translational  and  rotational  motion  described 
by  an  initial  angular  rotation  rate  (co),  initial  radial  velocity  (V,)  and  its  higher  order 
initial  angular  (y)  and  radial  acceleration  ( a ,)  and  so  on,  the  target  range  and  the 
rotation  angle  as  a  function  of  time  are  then 

R(t)  =  R0  +  Vtt  +—att2 +....  (2.2) 
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(2.3) 


1  9 

0{t )  =  0q  +  cot  +  —  yt  + _ 

Thus,  at  time  t,  (2. 1 )  becomes, 

Rp  =R(t)  +  x  cos  0(t)-y  sin  0(t)  (2.4) 

For  simplicity,  if  we  assume  the  reflectivity  function  of  target,  p(x,y),  to  be  a 
pure  function  of  the  target  geometry  and  that  it  does  not  change  with  the  above 
parameters,  then  the  base-band  of  the  returned  signal  from  the  scattering  point, 
P(x,y),  will  consist  simply  of  the  product  of  the  reflectivity  function  and  the  phase 
variation  due  to  the  distance,  RP\ 

Spit)  =  pix,  y)exp^-j2nfc  2R^t)^  =  p(x,  y)exy{-jyiRPt)}  (2.5) 


2.  Signal  Representation 

We  next  look  at  the  signal  representation  of  the  returned  echoes  for  the 
imaging  process.  At  any  given  time  t,  the  radar  receives  a  collection  of  returned 
signals  from  all  the  scattering  centers  of  the  target.  The  received  base-band 
signal  can  be  described  as  a  superposition  of  these  signals: 

-jin fc  I, tlxcly  (2.6) 

Substituting  (2.4)  into  (2.6),  we  obtain 

j^7if  °°  °° 

sR(t)  =  eJ  L  c  |  |  p(x,y)cxp{-j2n[xfx(t)  +  yfy(t)^dxdy  (2.7) 

-oo  -oo 

where 

/v(0  =  -^Lcos6*(r)  and  f  (t)  =  -^Lsm0{t) 
c  c 

Apart  from  the  pre-factor  (due  to  the  translational  motion)  in  (2.7),  the 
base-band  received  signal  at  a  given  time  can  be  seen  to  be  a  two-dimensional 
Fourier  transform  of  the  target  reflectivity  function.  This  is  the  basis  for  the 
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>(t)=  {  J  pix, y) ex p< 


method  of  extracting  the  target  reflectivity  information  by  applying  a  direct  two- 
dimensional  inverse  Fourier  transform  or  FFT  on  the  collected  data. 


3.  ISAR  Data  Collection  and  Processing 

In  almost  all  ISAR  systems,  the  signal  processing  is  performed  discretely. 
To  gain  more  intuition  into  this  procedure,  it  is  convenient  to  consider  the  discrete 
representation  of  the  base-band  signal  rather  than  its  continuous  form.  For  our 
earlier  discussion  on  ISAR  architecture,  we  set  the  frequency  dependence  for 
each  m  burst  of  pulses  to  increase  from  fc  in  steps  of  zl/from  pulse  to  pulse. 

fn  =  fc  +(n-l)Af  where  n  =  1,  2,  N  (2.8) 


Also,  the  base-band  signal  will  be  time-sampled  at  the  interval  corresponding  to 
the  pulse  interval  (T2)  plus  the  range  of  the  target  center.  Accordingly  (Jae, 
Thomas  and  Flores  [3]),  the  discrete  sampled  frequency  and  time  can  be 
represented  as 

2R 

t ni  n  =  \n  +  (m  ~  1)^V]T2  4 - -  where  m  =1,  2,  M  (2.9) 


Therefore,  (2.2),  (2.3)  and  (2.7)  can  be  expressed  in  terms  of  (2.8)  and 
(2.9)  as 


sR(m,n)  =  e  "  c  p(x,  y)  exp  {-y  2n  [xfx  (; m,n )  +  yfy  (m,  n)]}  dxdy 


2 R(m,n)  oo  oo 


R(m,n)  =  R0  +  V,t  nun  +  \attmJ  +.... 

2  f 

fM’n)  =  —cos(00  +  G)tmn+\ytmn-  +....) 


(2.10) 


2/„ 

c 


fv(m,ri)  = - -sm6(0o  +(otmn  +\yt,J- + 
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Figure  3  Row-column  decomposition  for  synthetic  ISAR  imaging. 


Figure  3  shows  the  image-processing  sequence  for  the  MxN  ISAR 
complex-valued  sampled  data  matrix  array.  The  image  construction  process 
begins  with  synthetic  range  profile  extraction  by  an  Inverse  Discrete  Fourier 
Transform  (IDFT).  This  operation  on  each  m  and  n  point  in  the  MxN  matrix  is 
given  as  such  (Jae,  Thomas  and  Flores  [3]). 


j  N -\ 

h  ( m  ,  n  )  =  —  Yj  sR(m,i)e 
N  i  =  0 


.2  n  ni 

1 - 

N 


where  n  =  1,  2,  N 


(2.11) 


The  resultant  MxN  matrix  is  analyzed  by  another  Fourier  process  -  the 
Discrete  Fourier  Transform  (DFT)  -  to  gather  the  cross-range  information.  If  we 
denote  the  output  of  the  transformed  point  in  the  resultant  MxN  matrix  by  D(m,n ), 
we  obtain 


J  M  - 1 

D{m,n )  =  —  V  h(k ,  n)e 
M  to 


.  2  n  km 
J~M 


where  m  =  1,  2,  M 


(2.12) 
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The  magnitude  of  the  elements  of  resultant  complex  matrix  forms  the  ISAR 
image. 

D.  ACHIEVING  HIGH  RANGE  RESOLUTION  THROUGH  STEPPED 

FREQUENCY  WAVEFORM  (SFW)  SIGNAL  PROCESSING 

Most  radars  today  attain  high  range  resolution  through  various  forms  of 
pulse  compression.  For  a  pulse  compression  radar,  the  achievable  range 
resolution  is  inversely  proportional  to  the  system  bandwidth,  and  thus,  to  gain 
high  resolution  the  entire  radar  system  —  from  transmitter  through  receiver  and 
including  data  collection  —  must  possess  the  large  bandwidth  associated  with 
the  desired  resolution.  Designing  a  radar  to  have  such  wide  bandwidth  at  all  its 
stages  can  be  an  expensive  and  a  rather  complex  issue. 

SFW  avoids  this  stringent  requirement  by  simply  varying  the  carrier 
frequency  prior  to  transmission  in  N  frequency  steps  associated  with  the  desired 
range  resolution.  The  same  stepped  carrier  frequency  is  also  used  as  the 
reference  to  the  mixer  to  down-convert  the  returned  signal  back  to  the  IF  stage. 
Consequently,  the  IF  stage,  the  detection  stage  and  later  processing  stages  no 
longer  require  large  bandwidth  to  maintain  the  resolution.  In  addition  to  relaxing 
the  wide  bandwidth  requirement  for  the  other  radar  system  stages,  the  SFW 
technique  also  has  the  advantage  that  the  returned  base-band  signal  is  a 
spectral  representation  of  the  target’s  reflectivity  with  range-invariant  magnitude. 
This  is  a  desirable  property  for  range  tracking  (Wehner  [4]). 

In  doing  so,  on  the  other  hand,  the  returned  signal  no  longer  preserves  its 
direct  phase  relationship  with  the  range  profile  of  the  target.  To  obtain  the  target 
range  profile  an  inverse  Fourier  transform  is  needed  to  get  back  the  target  range 
profile.  This  process  is  known  as  synthetic  range-profile  processing. 

1.  Synthetic  Range  Profile  Processing 

Figure  4  shows  an  SFW  transmission  and  its  associated  returned  echo 
pulse.  From  Wehner  [4],  each  received  echo  pulse  is  equivalent  to  a  replica  of 
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the  transmitted  SFW  with  a  time  delay  of  r(t)  =  2  (Rp  -  v,t)  /  c.  After  passing 
through  RF,  IF,  and  the  quadrature  detection  stage,  the  carrier  component  is 
removed  leaving  only  the  base-band  signal  with  a  phase  output  of  Wn(t)  =  2rfnT(t). 
This  signal  is,  in  turn,  sampled  at  S„  =  nT2  +  z>  +  2Rp/c  where  rr  is  the  receiving 
system  transfer  delay,  Rp  represents  the  range  of  a  point  target  (it  is  often 
estimated  through  range  alignment),  and  c  is  the  velocity  of  wave  propagation 
(Wehner  [4]). 


Transmitted  Pulses 


fc+(l-N)Af 


Figure  4  SFW  Pulse  and  echo  pulses. 
(After:  Figure  5.5,  pp  204,  Wehner  [4]) 


The  sampled  output  from  both  the  I-  and  Q-  channels  can  now  be  written 
in  complex  form  where  An  is  the  signal  amplitude  (which  is  usually  normalized  for 
ease  of  illustration  (Wehner  [4])): 

G„=Ane*»  (2.13) 


and 
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Wn  =  -2  nfn 


(2.14) 


2  R  2v 

f  2  R  V 

p  i 

nT2+rr  + 

c  c 

V  c  )_ 

The  N  sampled  output  data  for  each  burst  are  then  Fourier  transformed  by 
IDFT  or  FFT  to  obtain  a  series  of  synthetic  complex  range-profiles  of  the  target, 
Hi,  where  /  is  the  slant-range  position.  Following  Wehner  [4],  the  IDFT  is 
expressed  as 

N- 1 

0<l<N-\  (2.15) 

1=0 


Putting  (2.13)  and  (2.14)  into  (2.15),  and  assuming  zero  target  velocity,  the 
normalized  synthetic  response  of  the  target  is 


zj  V*  .  f  2;r  2  R 

H,=2^e*PJ  -T7h~2nf, — - 


;=0 


N 


(2.16) 


Replacing^  =fc  +  iAf  in  (2.16),  this  becomes 


H \  =  exp 


2  Rp  ^  ^ 

'ZexP  J 


J 


1=0 


2  ni 


N 


^  -2  NR  Af  ^ 
- p—  +  l 


(2.17) 


This  equation  can  be  further  simplified  using  the  same  analysis  as  for  antenna 
array  theory  (from  Wehner  [4])  to  become 


H ,  =  exp 


-j2nfc 


2  R. 


sin  7i  y 


2  sin 


ny 

N 


exp 


.N-l  2 nv 


\ 
V 

N~) 


(2.18) 


where 


~2NRpAf  t 
y  = - +1 

c 


(2.19) 


It  is  noted  for  (2.18)  that  the  resultant  response  represents  the  range  point 
spread  function  (PSF)  of  the  target  in  the  range  cell  index  I  where  the  highest 
peak  indicates  the  range  position  of  the  target  point.  The  collection  of  all  the 
range  cell  indices  from  I  =  0,  1,...,  N-l  is  known  to  be  the  synthetic  range  profile 
of  the  target. 
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2.  Unambiguous  Range  and  Range  Resolution 

From  both  (2.18)  and  (2.19),  the  peak  response  occurs  when  y  =  0,  +/-N, 
+/-2N,  +/-3N,...  corresponding  to  the  range  given  by  (Wehner  [4]) 

R  _  clo  c(lp  +  N)  c(/p  +  2N) , 
p  2NAf  ’  2NAf  ’  2NAf  K  ' 


The  unambiguous  range  length  is  then: 

A R  =  C{'- T  N) - — — or  c(/"  T  2N)  -  c(l°  T  N}  o,,..  =  -1-  (2.21) 

2N  Af  2NAf  2NAf  2NAf  2Af 

and  the  range  resolution  (using  Rayleigh’s  resolution  distance  criteria4) 

corresponds  to  the  width  of  the  PSF  of  each  target  point,  which  is 


AR 


P 


c 

2 m 


(2.22) 


Both  parameters  depend  on  the  frequency  step  size  but  the  range  resolution  can 
be  enhanced  by  placing  a  higher  number  of  pulses  within  each  burst. 


3.  Effect  of  Target  Velocity 

In  the  above,  the  synthetic  range  response  of  the  target  is  obtained  for  a 
target  with  zero  velocity.  However,  if  target  linear  velocity  is  present,  the 
synthetic  range  response  becomes 


2  R 


2v, 

f  ■  2RP  Y 

11 

t 

iT2  +  Tr  + 

l 

f 

c 

l  c  )_ 

JJ 

(2.23) 


In  this  case,  the  range  resolution  can  no  longer  be  simply  expressed  as  in  (2.22). 
Moreover,  analytically  deriving  the  range  resolution  ARP  from  (2.23)  is  also  not  a 
trivial  problem.  Through  numerical  solution,  Wehner  [4]  has  reportedly  shown 


4  Rayleigh  criterion  defines  the  resolution  between  two  point  spread  functions  to  be  the 
overlap  of  one  highest  maximum  peak  with  the  first  minimum  of  the  other  highest  peak.  For  a 
sine  function,  this  occurs  at  the  first  null  in  which  the  argument  in  the  sine  function  is  n. 
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that  the  higher  velocity  can  cause  the  width  of  the  PSF  to  broaden,  thus  resulting 
in  poorer  range  resolution. 


4.  Range  Offset  and  Range  Walk 

The  presence  of  nonzero  translational  target  velocity  has  the 
consequence  of  creating  range  offset  and  range  walk  effects,  and  therefore 
distorts  the  radar  image.  When  ISAR  constructs  a  radar  image  from  the  collected 
data,  it  presumes  that  the  range  of  every  target  scattering  point  does  not  vary 
during  the  integration  time  interval.  Range  offset  and  range  walk  result  when  this 
presumed  condition  fails  to  be  true  due  to  the  presence  of  radial  velocity 
components.  The  extent  of  the  range  walk  and  range  offset  can  be  quantified  by 
the  number  of  cells  (AN)  that  the  range  is  shifted  over  the  entire  time  integration 
period  ( T)  for  the  radar  time  delay  resolution  of  At  =  1  / NAf.  Thus,  from  Wehner 
[4],  AN  is  expressed  as 

AN  =  —  =  NA/St  (2.24) 

At 


where  St  represents  the  time  delay  shift  of  the  target  due  to  the  presence  of 
radial  velocity  over  the  period  where  St  is  evaluated.  This  delay  shift  is 
determined  from  the  total  group  delay  and  can  be  derived  by  taking  the  derivative 
of  (2.14)  and  subtracting  the  time  delay  2Rp/c.  In  general,  Aris  cumulative  over 
the  number  of  bursts,  m,  and  is  given  by  Wehner  [4]  as 


St  = 


2v, 


cAf 


( 


fcT2+  A/ 


27? 

r.  H - b  NmT-, 


(2.25) 


Putting  (2.25)  into  (2.24),  and  noting  that  NmT2  »  t ,  +  2R  /  c  is  true  in  most 
cases,  we  obtain 

AW  ~  -  (fc  +  NmAf )  vt  (2.26) 

For  the  entire  image  frame  time,  we  let  m  =  M,  and  obtain 
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(2.27) 


A N~- 


2 (ivr2)  |  2M (NT2)N Af 


The  first  term  in  (2.27)  represents  the  range  offset  which  increases  with 
higher  values  of  N,  higher  carrier  frequency  and  lower  PRF.  The  second  term  is 
the  result  known  as  range  walk  and  is  cumulative  from  profile  to  profile  (burst  to 
burst)  depending  on  size  of  M.  Unlike  range  offset,  the  non-constant  range  shift 
due  to  the  cumulative  nature  of  range  walk  tends  to  move  the  same  scattering 
points  out  of  the  range-cell  column  (known  as  range  cell  migration)  and  affects 
the  cross-range  information  extraction  process  (which,  in  effect,  takes  DFTs  of 
each  range-cell  column).  This  will  generally  result  in  distortion  of  the  ISAR  image. 
The  choice  of  Af  also  contributes  to  the  range  walk  effect.  For  better  range 
resolution,  we  require  high  Af  but  range  walk  will  tend  to  be  amplified  by  higher 

4f- 

To  correct  the  distortion  caused  by  range  cell  migration  of  the  scattering 
points,  the  phase  error  due  to  the  target  velocity  has  to  be  mitigated.  The 
simplest  method  for  eliminating  this  error  is  to  multiply  the  sampled  output  data 
by  an  equivalent  complex  factor  (Wehner  [4]) 


g  =  exp 


-fin  ft 


2v. 


f  ij  5 

M  iT2+Tr+-?- 

V  '  c 


(2.28) 


where  vt  and  R  are  estimated  from  the  Doppler  and  range  tracking  of  the 

sampled  output  data.  This  is  known  as  translational  motion  compensation. 
However,  it  is  evident  from  (2.28)  that  the  effectiveness  of  this  kind  of 
compensation  relies  heavily  on  the  precision  of  Doppler  and  range  tracking 
techniques  in  order  to  get  an  accurate  estimated  of  the  target  range  and  its  radial 
velocity.  Compensation  methods  get  more  complex  if  the  target  velocity  varies 
significantly  within  the  integration  period. 

There  are  many  methods  today  that  can  yield  accurate  motion  estimation. 
One  of  the  known  methods  is  the  entropy  measurement  but  this  method  is 
computationally  intensive  (Jae  and  Thomas  and  Flores  [3]).  Other  techniques 
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include  the  Cramer-Rao  Burst  Derivatives  Approach  (Jae  and  Thomas  and 
Flores  [3]),  the  Phase  Difference  method  (Wehner  [4]),  the  least  square  motion 
parameter  estimation  (Jae  and  Thomas  and  Flores  [3]),  and  even  time-frequency 
analysis  approaches  (Jae  and  Thomas  and  Flores  [3]  and  Chen  [2]). 


E.  ACHIEVING  HIGH  CROSS-RANGE  RESOLUTION  USING  SYNTHETIC 

APERTURE 

1 .  ISAR  Theory  from  an  Aperture  Viewpoint 

The  method  by  which  ISAR  achieves  high  cross-range  resolution  is  not 
much  different  than  that  of  conventional  SAR  and,  in  both  cases,  the  bounds  are 
linked  closely  to  the  angular  extent  over  which  the  radar  views  the  target.  For  an 
unfocussed  SAR,  it  is  the  maximum  linear  path  length  before  the  quadratic  phase 
error  becomes  significant  to  create  the  defocusing  effect  and  degradation  of 
resolution.  (This  phase  error  results  from  the  difference  between  the  linear  path 
and  the  curvature  of  the  range  to  a  point  target.)  In  the  case  of  ISAR,  the  same 
behavior  holds  except  that  the  path  length  is  now  a  function  of  the  target  angular 
rotation  rate  and  the  dwell  time  (which  is  also  usually  the  image  coherent 
processing  or  integration  time). 

Consider  a  stationary  radar  viewing  a  rotating  target  (see  Figure  5(a))  over 
an  angle  segment  W  (aperture)  over  the  integration  period  (we  assume  that  the 
translational  motion  is  perfectly  compensated).  This  geometry  can  be  viewed  as 
equivalent  to  the  radar  (like  SAR)  in  motion  observing  a  non-rotating  target  as  it 
passed  the  same  angle  segment  as  shown  in  Figure  5(b).  We  assume  that  the 
radar  is  tracking  the  target  in  the  far-field  so  that  the  angular  extent  of  the 
aperture,  projected  as  path  length,  L,  is  L  «  R.  Also,  the  length  of  the  observed 
target,  v,  is  generally  y  «  L.  In  such  cases,  from  Figure  5(c),  the  two-way  phase 
advance  of  the  target  scattering  point,  <p,  at  azimuth  position,  from  the  radar 
boresight  is 

Air 

cp{x)~2kSr  = — xsin  <f>y  (2.29) 

X 
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Stationary  Radar 


(a)  ISAR  (After:  Figure  7.3(a),  pp  344,  Wehner  [4]) 


(b)  SAR  equivalent  (After:  Figure  7.3(b),  pp  344,  Wehner  [4]) 


(c)  Approximation  for  L  «  R  and  y  «  L  (After:  Figure  7.3(c),  pp  344,  Wehner  [4]) 

Figure  5  SAR  and  ISAR  equivalence. 
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The  term  sin  fa  approaches  y  /  R  for  y  «  R,  and  the  projected  length  x 
approaches  Rtfr  under  small  angle  approximation.  Hence,  (2.29)  can  be 
expressed  in  terms  of  ^and  becomes  <p(<f)  =  4n  (fry  /  X.  Using  similar  arguments 
as  in  linear  antenna  array  theory,  the  normalized  integrated  response  over  the 
integration  angle,  lF,  behaves  like  a  point  spread  function  (PSF)  about  the  cross¬ 
range  position,  v: 

Z(y )  =  \\  exp (j^(fry)d(fr  =  Sinf  /:  ¥y^  y  (2.30) 

j-t  if  wy 


Using  the  Rayleigh  resolution  criteria,  we  can  determine  the  cross-range 
resolution  between  two  point  targets  to  be  Arc  =  X  /  2W.  For  small  integration 
angle  ISAR,  we  set  {F=  coT,  where  co  denotes  the  target  angular  rotation  rate  and 
Tis  the  integration  time.  Then,  the  cross-range  resolution  is  given  by 


A  r  =  ■ 


X 


2  coT 


(2.31) 


Like  SAR,  the  size  of  the  target  imposes  a  maximum  allowable  W  for 
which  this  cross-range  resolution  can  be  achieved.  It  can  be  shown  (see  Wehner 
[4])  by  setting  the  maximum  allowable  phase  derivation  to  be  no  more  than  tc /  8 
rad  as  the  criteria  for  maintaining  focus,  that  the  maximum  integration  angle 
before  defocusing  occurs  is 


V 


max 


(2.32) 


where  the  maximum  target  length  is  L,  and  to  maintain  focus  it  is  necessary  that 


L< 


{Arcf 

2X 


(2.33) 


2.  Range  Doppler  Imaging  and  Cross  Resolution 

The  cross-range  resolution  in  the  above  situation  is  obtained  by  drawing 
an  analogy  to  SAR  in  terms  of  geometry.  The  aim  is  to  emphasize  the 
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importance  of  maintaining  a  small  integration  angle,  %  within  the  linear  regime 
for  which  the  quadratic  phase  is  not  significant.  Unlike  SAR,  however,  in  ISAR 
this  angle  is  affected  by  the  target  rotation  rate  and  is  often  not  a  known  priori  — 
nor  can  it  be  generally  measured  directly  by  the  radar. 

In  most  —  if  not  all  —  ISAR  systems,  Doppler  processing  techniques  are 
directly  applied  to  the  signal  to  obtain  the  cross-range  information.  This 
technique  avoids  the  need  for  exact  knowledge  of  the  target  rotation  rate,  co,  as 
long  as  the  product  of  coT  remains  within  the  small  integration  angle  limit.  The 
relationship  between  target  rotation,  the  point  scatterer’s  cross-range  position 
and  its  resulting  Doppler  shift  can  be  understood  by  referring  to  Figure  6.  If  we 
ignore  the  Doppler  component  due  to  translational  motion  (it  will  be  shown  later 
that  this  component  is  just  a  linearly  added  term  to  the  Doppler  shift),  then  the 
Doppler  shift  due  to  the  instantaneous  velocity  of  the  scattering  point  at  a 
distance  rc  from  the  center  of  target,  arc,  is 


Figure  6  Radial  velocity  from  a  point  scatterer  on  a  rotating  target. 
(After:  Figure  7.6,  pp  351,  Wehner  [4]) 


From  Equation  (2.34)  (after  some  manipulation),  the  cross-range  resolution  is 
seen  to  be  directly  related  to  the  Doppler  resolution  of  the  radar,  AfD. 
Furthermore,  since  AfD  =  1/T,  the  cross-range  resolution  can  be  expressed  as 


A  r 


/L  A 

- 4 fn  = - 

2(0  2  coT 


(2.35) 


24 


This  expression  is  the  same  as  that  obtained  earlier  for  the  cross-range 
resolution. 

3.  Doppler  Frequency  Shift  Due  to  Target  Motion 

The  simple  relationships  above  may  not  reveal  much  about  the  effects  of 
target  motion  of  the  ISAR  image.  To  give  further  insight,  let’s  consider  the  phase 
relationship  derived  earlier  in  Equation  (2.5).  An  approximation  of  the  Doppler 
shift  induced  by  the  target’s  motion  can  be  determined  by  taking  the  time- 
derivative  of  the  phase  term  y/(Rpt)  where  Rpt  =R(t)  +  x  cos  0(t)-y  sin  9(t) , 

R(t)  =  R0  +  Vtt  +  ^a,t2  +....  and  0(t)  =  90  +cot+\yt2  + .  Neglecting  the  second 

order  and  higher  terms  of  R(t)  and  9(t),  we  obtain 
d  2  f  2  f 

f D  = — i//{RPt)  =  ^-£-Vt  +  {-x&>[sin  9(]  cos  cot  +  cos  9f)  sin  coti 

dt  c  c 

-yco\cos90  cos  cot  +  sin  0O  sin<xv]}  (2.36) 

f  DTrans  foRot 

The  overall  Doppler  shift  can  be  seen  as  a  summation  of  the  Doppler  shift 
due  to  translational  motion,  form™  =  2fcVt/c,  and  that  due  to  rotational  movement, 
fDRot,  determined  by  the  remaining  terms.  We  can  further  expand  sin  cot  =  cot  -  co3t 3 
/  6  +  ...  and  cos  cot  =  1  -  art2  /  2  +  ...  and  ignore  the  higher  order  terms  since  the 
product  of  both  co  and  t  is  generally  very  small,  cot3  «  1  and  cot2  «  1.  Thus, 
foRo,  becomes 

foRot  =—^L{-o)[xsin90+ycos90\-o1t[xcos90-ysin90§  (2.37) 

The  first  term  in  (2.37)  is  used  for  ISAR  imaging.  But  there  exist  quadratic 
phase  terms  even  though  the  rotation  rate  remains  constant.  These  quadratic 
terms  vary  with  time  and  serve  to  introduce  an  additional  phase  error  that  can 
defocus  and  distort  the  image  during  its  construction  process. 


25 


4.  Distortion  Produced  by  Target  Rotation 

Unlike  errors  due  to  translational  motion,  the  phase  error  due  to  target 
rotation  can  cause  cell  migration  in  both  range  and  cross-range  dimensions  and, 
hence,  distortion  in  the  ISAR  image.  From  the  2D  radar-target  geometry  given 
earlier,  cell  migration  in  the  range  dimension  due  to  target  rotation  is  most 
significant  when  the  initial  target  angle  is  90  =  it/ 2.  The  number  of  cells  moved  in 
the  range  dimension  over  the  integration  period  T  is  then  (following  Wehner  [4]) 

Rp  ( t )  =  R0  +Vtt  +  x  cos(6*0  +  cot)  -  y  sin(6*0  +  cot) 


To  account  for  time  delay  shift  due  to  rotation  alone,  we  set  Vt  =  0. 
Expanding  the  sine  and  cosine  terms  in  RP  and  setting  Q0  =  n  /  2,  we  can 
approximate  AM’  as 


AM'  = 


1 

A  r 

J_ 

A  r 


{y-x  sin  coT  -  y  cos  coT } 


\y-x 

-y 

1 

6 

L  2  JJ 

xcoT 
A  r 


(2.39) 


Cross-range  cell  migration  is  most  significant  when  90  =  0.  The  number  of 
cross-range  cells  by  which  the  scattering  point  is  offset  from  its  correct  position 
due  to  target  rotation  rate  over  the  integration  period  T  is  then 

AM  =abs^-{fDRm((,)-fmJT)}  (2.40) 

A Td 


Substituting  (2.37)  into  (2.40)  and  letting  90=  0  and  AfD  =  1  /  T,  we  get 


AM  = 


2  UcoTfx 


(2.41) 


Earlier,  we  defined  the  cross-range  resolution  as  Arc  =  A./ 2coT  =  c / 2fccoT,  and  so 
by  substituting  this  result  into  the  above  equation,  we  obtain  the  same  expression 
as  in  Equation  (2.39). 
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(2.42) 


A M  = 


xcoT 
A  r 

c 


The  extent  of  the  cell  migration  is  dependent  on  the  target  size,  its  rotation 
rate  and  the  integration  time.  Also,  higher  resolution  requirements  over  the  same 
period  T  would  mean  increasing  the  number  of  range  (N)  and  cross-range  cell 
(M),  and  a  corresponding  increase  in  AM  and  AM’.  This  is  similar  to  the  effect  of 
increasing  the  range  resolution  discussed  in  the  previous  section. 

Often,  ISAR  can  be  made  to  remain  focused  over  the  entire  rotation  angle 
W  if  the  target  phase  drift  in  both  range  and  cross-range  (also  known  as  blur 
radius)  caused  by  target  rotation  is  constrained  to  lie  within  one  cell.  From 
Equation  (2.42),  if  we  set  AM  <  1  for  focused  ISAR  and  for  coT  =  X  /  2Arc,  this 
yields 


x< 


2(A  ;;)2 
X 


(2.43) 


Equation  (2.43)  is  almost  the  same  as  Equation  (2.33)  which  was  obtained 
by  setting  the  maximal  phase  deviation  to  n/ 8.  The  equation  is  determined  by 
“the  extremes  of  the  target  extent”  (Wehner  [4]),  the  required  resolution  and  the 
frequency  used. 
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III.  BLAST  LOADING  ON  AIRCRAFT 


A.  BLAST  IMPACT  ON  AIRBORNE  TARGETS 

Explosive  munitions,  guided  or  unguided,  are  extensively  used  in  most 
weapons  today  against  all  kinds  of  targets:  aircrafts,  ships,  tanks,  troops, 
infrastructures,  etc.  These  munitions  are  attractive  because  of  the  enormous 
chemical  energy  that  is  released  in  a  very  short  period  when  the  explosive  is 
detonated.  The  huge  energy  deposited  into  the  medium  creates  a  high  velocity 
high  pressure  shock  front  that  has  damaging  effects  on  the  target  structure  — 
especially  for  buildings  and  fixed  infrastructures.  When  used  against  airborne 
targets,  the  primary  goal  is  to  make  use  of  this  high  velocity  pressure  wave  to 
launch  small  pre-fa bricated  fragments  and  shrapnel  which  will  penetrate  different 
parts  of  the  target  body,  causing  a  variety  of  damage  mechanisms  such  as 
structural  tearing,  explosion  of  the  fuel  tank,  loss  of  flight  control  and  electronics, 
and  injuring  or  killing  the  pilot. 

Owning  to  its  finite  inertial  mass,  the  speed  at  which  these  fragments  and 
shrapnel  travel  will  in  general  be  less  than  the  shock-wave  speed.  As  a 
consequence,  the  blast  wave  will  interact  with  the  target  body  before  the 
fragments  and  shrapnel  arrive.  Most  modern  airborne  targets  are  constructed 
from  material  with  tensile  strength  able  to  sustain  high  aerodynamic  and 
environmental  drags  of  flight.  Also,  because  of  the  degrees  of  freedom  of  the 
aircraft  in  air,  the  net  effect  of  a  blast  force  is  more  likely  to  cause  translational 
and  rotational  motions  (as  if  it  is  treated  as  a  rigid  body)  than  yielding  its  structure 
strength.  A  brief  impulse  moment  can  thus  be  captured  by  ISAR  for  Battle 
Damage  Assessment  (BDA)  purposes. 

For  the  fragments  and  shrapnel  that  penetrate  the  target  body,  different 
damage  mechanisms  can  occur  depending  on  the  nature  of  the  penetration.  If 
the  fragments  hit  a  fuel  tank,  a  secondary  explosion  will  occur  and  may  cause 
further  structural  breakup.  Alternatively,  the  fragments  may  cause  the  target  to 
lose  flight  control  by  damaging  the  electronics,  aero-structures,  control 
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hydraulics,  or  even  killing  the  pilot.  In  such  cases,  the  target  will  display  erratic 
maneuvers.  Other  scenarios  include  tearing  out  the  tail  and/or  the  wings, 
damaging  the  engine,  detonating  the  weapons  stores,  etc.  Although  these 
effects  are  also  captured  by  ISAR,  they  are  in  general  difficult  to  completely 
account  for  even  in  computer  modeling  and  simulations. 

In  this  thesis,  we  would  like  to  briefly  examine  the  effects  of  an  explosive 
blast  on  the  motions  of  the  airborne  target  and  demonstrate  how  it  can  be 
detected  using  ISAR.  These  imposed  motions  are  highly  dynamic  and  depend 
on  strength  of  the  explosion,  the  relative  distance  and  direction  in  which  the 
detonation  occurs,  and  the  surface  area  of  the  aircraft  subjected  to  the  blast 
loading.  Very  often,  the  solution  to  these  complex  motions  can  not  be  easily 
derived  analytically.  However,  numerical  solution  of  these  complex  equations  of 
motion  is  possible.  Hence,  this  chapter  aims  to  provide  discussions  linking  the 
theory  of  explosive  blast  in  air,  the  loading  of  the  blast  force  on  the  aircraft,  and 
the  aircraft’s  translational  and  rotational  responses  through  a  set  of  complex 
equations  of  motion.  In  addition,  a  methodological  approach  is  also  proposed  to 
create  a  computer  model  to  solve  these  complex  equations  numerically. 

B.  CHARACTERISTIC  OF  EXPLOSIVE  BLAST  IN  AIR 

1.  Formation  of  the  Blast  Wave 

Blast  waves  are  formed  as  a  result  of  the  ambient  atmosphere  being 
forcibly  pushed  back  by  the  high  pressure  gases  produced  from  an  explosion.  In 
the  initial  formation,  the  shape  of  the  pressure  pulse  from  a  conventional 
chemical  explosion  is  highly  dependent  on  the  geometric  construction  of  the 
explosive  charge.  However,  as  the  pulses  travel  in  the  medium,  the  higher 
pressure  portions  will  possess  higher  speed  than  the  lower  pressure  parts 
(Kinney  and  Graham  [10]).  Progressively,  the  wave  front  becomes  steeper  as 
the  pressure  pulse  move  out  and  develops  a  discontinuity  known  as  an  explosive 
shock  as  illustrated  in  Figure  7.  Under  this  condition,  increasing  amounts  of 
energy  will  be  needed  to  further  alter  the  shock  front  shape.  The  front  will  largely 


30 


remain  in  this  form  and  acts  as  an  advancing  blast  wave  system  propagating 
away  from  the  charge  in  the  radial  direction. 


(a)  (b)  (c) 

Figure  7  Development  of  explosive  shock,  (a)  Initial  pressure  pulse,  (b)  successive 
configurations,  and  (c)  Final  form  when  fully  developed. 

(After:  Fig  6-1 ,  pp  89,  Kinney  and  Graham  [10]) 


It  is  interesting  to  note  that  “the  fully  developed  explosive  blast-wave 
system  is  always  formed  with  about  the  same  general  configuration  no  matter 
what  is  assumed  for  the  initial  positive  pulse.  That  is,  any  initial  configuration  is 
lost  and  all  blast  waves  at  reasonable  distances  from  the  center  of  an  explosion 
show  similar  wave  forms”  (Kinney  and  Graham  [10]). 

The  amount  of  energy  released  by  the  chemical  reaction  process  from  the 
explosion  is  an  important  factor  in  determine  the  intensity  of  the  initial  shock  front 
and  its  subsequent  impulse  shape.  In  its  initial  form,  the  explosion  process 
entails  a  vanishingly  small  volume  of  gases  under  high  pressure.  As  the  gases 
expand  outward,  the  same  amount  of  energy  will  need  to  overcome  a  larger 
volume  of  the  atmosphere  (which  acts  to  resist  the  expansion  process)  hence 
lowering  its  resultant  shock  front  intensity  and  ‘stretching’  the  wave  shape  as 
shown  in  Figure  8. 


31 


Figure  8  Typical  pressure-distance  curves  for  successive  times  after  an  explosion. 

(After:  Fig  6-2,  pp  90,  Kinney  and  Graham  [10]) 


At  some  appreciable  distance,  the  shape  of  the  blast  wave  will  depict  a 
series  of  positive  and  negative  phases  of  the  pressure.  In  general,  the  first 
positive  pressure  phase  is  far  more  intense  in  causing  damage  than  the  negative 
ones  and  its  subsequent  rarefactions  which  are  limited  in  magnitude  to  those  of 
the  atmosphere  (Kinney  and  Graham  [10]). 


2.  Characteristic  of  the  Blast  Wave 


Figure  9  Typical  pressure-time  curves  for  explosive  blast  wave. 
(After:  Fig  6-5,  pp  99,  Kinney  and  Graham  [10]) 
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Figure  9  shows  a  typical  pressure-time  history  characteristic  curve  of  an 
explosive  blast  at  a  certain  location  away  from  the  center  of  the  explosion  (after  it 
is  fully  developed).  At  arrival  time  ( ta )  after  the  explosion,  there  exists  a  sudden 
jump  in  the  pressure  intensity  (known  as  peak  overpressure)  from  its 
atmospheric  level.  An  object  at  this  location  will  be  subjected  to  an 
instantaneous  lateral  force  equal  to  the  product  of  the  overpressure  and  the 
projected  area  in  the  plane  of  the  blast  wave.  However,  this  is  not  a  stable 
condition,  and  the  overpressure  at  the  location  will  begin  to  decay  exponentially 
as  the  volume  of  gases  continues  to  expand  outward,  leading  to  smaller  negative 
and  subsequently  smaller  positive  rarefaction  phases.  For  most  purposes,  the 
main  interest  lies  in  the  first  positive  phase  of  the  blast  wave  where  the  damaging 
effect  is  considerably  more  significant  than  its  later  ones. 

There  are  four  independent  parameters  that  can  be  used  to  completely 
describe  the  first  positive  phase  of  blast  wave  shape:  (a)  the  initial  shock  intensity 
specified  by  peak  overpressure  ( pa ),  (b)  the  duration  of  the  positive  phase  of  the 
blast  (td),  (c)  the  impulse  (force-time  product)  per  unit  area  (I/A),  and  (d)  the 
arrival  time  (ta)  which  merely  introduces  an  overall  time  shift  of  the  pressure-time 
history  curve. 

An  extension  of  the  logarithmic  decay  expression  from  Kinney  and 
Graham  [10]  linking  these  parameters  to  the  blast  wave  shape  can  be  given  as 


p(t)  = 


Po 

0 


f  t-t  A 
V  fd  J 


exp 


a(t-ta) 


where  ta  <t  <td 
otherwise 


(3.1) 


where  p  is  the  instantaneous  overpressure  at  time  t,  and  a  is  the  waveform 
parameter  related  to  the  impulse  per  unit  area. 
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3.  Peak  Overpressure 

Theoretical  analysis  of  the  peak  shock  overpressure  uses  the  same 
approach  as  for  a  normal  shock.  It  is,  however,  complicated  by  the  spherical 
divergence  of  the  shock  front  and  its  transient  nature.  Both  theoretical  studies  (in 
some  cases,  through  the  use  of  computer  solutions)  and  experimental 
measurements  are  usually  correlated  to  obtain  a  reasonable  expression  relating 
peak  intensity  to  the  distance.  Following  Kinney  and  Graham  [10],  the 
empirically  derived  peak  overpressure-to-distance  relationship  can  be  expressed 
as  the  ratio  of  the  overpressure  to  the  standard  atmospheric  pressure  of 
1013.25mbars  at  15deg  C. 
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(3.2) 


where  p,/Pa  is  overpressure  ratio  and  Z  is  the  scaled  distance  of  a  charge  weight 
at  lkg  TNT  equivalent.  The  scaling  of  Z  to  the  actual  distance  will  be  discussed 
below. 


4.  Arrival  Time 

From  Kinney  and  Graham  [10],  we  noted  that  the  velocity  of  the  shock 
front  is  uniquely  related  to  po/Pa  which  in  turn  depends  on  the  distance  between 
the  observed  point  and  the  center  of  the  explosion.  Hence,  given  the  peak 
overpressure  and  the  distance  over  which  it  occurs,  we  can  obtain  the  arrival 
time  by  noting 


u 


x 


dr 

dt 


=  aM„ 


(3.3) 


where  ux  is  the  shock  front  velocity,  Mx  is  its  Mach  number,  r  is  the  radial 
distance,  and  ax  is  the  speed  of  sound  in  the  undisturbed  atmosphere.  Solving 
(3.3)  by  separation  of  variables  and  integrating  the  Mach  number  from  the  center 
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of  the  charge  (assumed  to  be  rc  =  0)  to  the  distance  the  shock  front  arrived,  we 
obtain 


t„  = 


=  —  \—dr 

n  J  M. 


~x  0 


(3.4) 


where  the  expression  for  shock  velocity  in  term  of  Mach  number  (given  by  Kinney 
and  Graham  [10])  is 

M  =  1  +  — +  where  k  =  1.4  (3-5) 

2  kPa 


5.  Blast  Duration 

The  extent  of  the  damage  caused  by  the  blast  wave  is  also  closely  linked 
to  the  duration  over  which  the  force  is  applied.  Since  the  impulse  (equivalent  to 
the  area  under  the  blast  wave  shape)  of  the  positive  phase  is  much  higher  than 
its  subsequent  phases,  it  is  often  used  as  an  index  for  the  blast  duration  although 
in  “some  cases  the  negative  phase  duration  can  be  twice  as  long”  (Kinney  and 
Graham  [10]).  Thus,  the  blast  duration  at  an  observed  point  is  defined  to  be  “the 
time  between  the  passing  of  the  shock  front  and  the  end  of  the  positive  pressure 
phase  as  marked  by  a  zero  overpressure”  (Kinney  and  Graham  [10]). 

The  zero  overpressure  condition  is  a  characteristic  of  the  sound  wave 
since  the  shock  front  velocity  will  decrease  to  the  speed  of  sound  in  the  medium. 
Hence  the  blast  duration  is  also  dependent  on  the  shock  velocity.  Using  a  similar 
approach  to  that  for  deriving  the  arrival  time,  Kinney  and  Graham  [10]  give  the 
relationship  between  the  blast  duration  and  the  distance  to  which  the  shock  front 
begins  to  develop  as 
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where  the  term  t/  WI/3  is  the  duration  in  msec  for  a  1  kg  TNT  explosion  and  Z  is 
the  same  scaled  distance  as  in  the  peak  overpressure  expression.  The  actual 
blast  duration  can  be  obtained  by  multiplying  (3.6)  by  W1/3  where  W\s  the  weight 
of  the  explosive  used  in  kg. 


6.  Impulse  Per  Unit  Area 

Impulse  has  “the  dimensions  of  the  force-time  product”  and  can  be 
obtained  graphically  from  the  area  under  the  first  positive  phase  of  the  blast  wave 
characteristic.  This  factor  determines  the  extent  of  damage  to  the  target 
structure.  From  Figure  9,  we  can  see  that  —  apart  from  the  peak  overpressure 
and  blast  duration  —  the  rate  of  decay  of  the  overpressure  (known  as  the 
waveform  parameter,  a)  also  affects  the  amount  of  impulse  force  acting  on  the 
target.  A  rapid  decay  blast  wave  with  the  same  intensity  will  have  less  damaging 
effect  on  the  target  than  the  slow  decaying  ones. 


However,  it  is  often  difficult  to  determine  a  from  analytical  analysis  or 
direct  measurement.  In  contrast,  the  impulse  per  unit  area  (I/A)  can  be  easily 
obtained  by  experiment,  and  the  empirical  relationship  between  I/A  and  the 
scaled  distance  (from  Kinney  and  Graham  [10])  can  be  expressed  as 
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(3.7) 


With  Equation  (3.7),  the  waveform  parameter  in  turn  can  be  determined  by 
rearranging  the  following  integral  expression 
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7.  Hopkinson  Scaling  Law 

Scaling  of  parameters  allows  the  quantitative  characteristic  of  any 
explosion  (especially  those  drawn  from  experimental  data)  to  be  used  for 
solutions  to  more  general  explosive  blast  wave  problems.  The  scaling  law  is 
fundamentally  based  on  the  geometrical  similarity  between  the  reference  and 
actual  object.  The  general  principle  for  explosive  scaling  is  to  consider  the  two 
objects  to  be  spherical.  Then,  by  geometry,  the  ratio  of  their  volumes  is 
proportional  to  the  third  power  of  the  ratio  of  their  diameters,  and  if  the  two 
objects  have  uniform  density  distribution,  their  mass  ratio  also  follows  the  same 
third  power  rule. 

By  extending  the  above  argument,  the  blast  wave  of  two  explosions  can 
be  said  to  be  identical  if  the  ratio  of  their  distance  is  proportional  to  the  cube  root 
of  their  energy  release.  If  we  take  the  medium  into  account,  it  has  been  shown 
by  Kinney  and  Graham  [10]  that  the  scaling  for  distance  can  be  expressed  as 
follows: 


(< actual  _  dis  tan  ce){atmo spheric  _  density ) 


1/3 


(< energy  _  release) 


1/3 


=  ( scaled  _  dis  tan  ce) 


(3.9) 


Since  the  energy  released  is  related  to  the  velocity  of  the  shock  front, 
which  in  turn  depends  on  the  peak  overpressure  ratio  and  hence  the  atmospheric 
density,  the  above  expression  can  be  simplified  to  become 


Z  =  ( scaled  _  dis  tan  ce)  = 


fd  x  ( actual  _  dis  tan  ce) 
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1/3 


(3.10) 


where  fd  is  known  as  the  atmospheric  transmission  factor  and  accounts  for  the 
difference  between  the  actual  and  standard  reference  atmosphere  (designated 
by  subscript  o).  fd  can  be  found  using 
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The  overpressure  is  indirectly  determined  by  the  scaled  distance  relation, 
and  the  actual  overpressure  is  simply 

Overpressure  =  overpressure  _  ratio  x  atmospheric  _  pressure  (3.12) 

This  is  the  situation  of  the  Hopkinson  Scaling  Law  (Ref  from  Kinney  and  Graham 

[10]). 

The  scaling  law  for  time  can  also  be  derived  from  the  scaled  to  actual 
distance  ratio  by  noting  that  the  time  associated  with  the  blast  is  the  integration  of 
the  speed  of  the  blast  wave  over  the  distance  traveled.  Thus,  by  Kinney  and 
Graham  [10],  the  actual  time  is  related  to  the  scaled  time  by 

,  ,  W113 (scaled  time )  .  0. 

(i actual  _tune)  = - = -  (o .  1  3 ) 

ft 


and  ft  is  the  transmission  factor  for  time  taking  also  into  account  the  ratio 
between  the  shock  front  speed,  a,  and  the  sound  speed  in  the  medium,  a0.  The 
expression  for  f  is  given  as: 
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C.  DYNAMIC  LOADING  OF  BLAST  WAVE  ON  AIRCRAFT 

1.  Translational  and  Rotational  Responses 

An  object  placed  in  the  field  of  an  explosive  blast  wave  will  experience  a 
dynamic  load  characterized  by  the  sudden  increase  in  its  peak  pressure  value 
followed  by  a  gradual  decaying  process.  The  net  effect  depends  on  the 
distribution  of  the  pressure  wave  across  the  aircraft’s  surface,  which  depends  on 
the  orientation,  geometrical  shape  and  construction  of  the  object,  and  its  resistive 
forces  such  as  bending  moments,  drag  and  aircraft’s  thrust  apart  from  its  inertia. 
The  end  result  is  most  likely  a  non-symmetrical  distribution  of  the  pressure 
intensity  across  the  exposed  area  and  this  loading  also  changes  with  time.  Thus, 
it  is  generally  a  difficult  problem  to  solve  analytically. 
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In  this  thesis,  the  main  interest  is  to  examine  the  translational  and 
rotational  responses  of  the  aircraft  as  a  whole.  To  accomplish  this  goal  we 
assume  the  aircraft  to  be  a  rigid  body  with  center  of  mass  located  at  the  origin  of 
the  aircraft  coordinates  (as  in  Chapter  II  except  that  we  will  expand  the  case  into 
the  third  dimension).  Consider  the  geometry  of  the  three  dimensional  (3D) 
problem  as  shown  in  Figure  10.  We  assume  the  blast  wave  originates  from  a 
point  source  and  diverges  outward  to  the  target.  Different  parts  of  the  aircraft 
surface  area  exposed  to  the  blast  wave  will  then  experience  different  peak 
overpressure  intensities  of  different  durations  and  different  arrival  times  (because 
of  their  relative  distances  to  the  blast  point). 


x 


Figure  1 0  3D  geometry  of  the  interaction  of  blast  wave  on  the  aircraft  surface. 

A  directed  area  element  on  the  aircraft  surface  at  point  P  can  be  described 
by  the  usual  area  element  (AA)  and  normal  vector,  «„  describing  the  orientation 
of  the  aircraft  surface.  The  point  P  is  denoted  by  the  vector  a ,.  A  spherical 
diverging  blast  wave  will  travel  in  the  radial  direction  along  the  directed  line 
linking  the  blast  point  to  the  point  where  the  pressure  interacts  with  the  target 
(point  P).  Denote  this  vector  by  RBi,  which  can  be  derived  if  we  know  a„  and  the 
distance  between  the  center  of  mass  of  the  aircraft  and  the  blast  point,  denoted 
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by  the  vector,  Rc.  This  distance  can  be  taken  to  be  the  miss  distance  of  the 
target  (if  we  assume  that  the  proximity  fuse  of  the  explosive  munitions  triggers  at 
the  point  where  the  largest  surface  area  of  the  target  is  sensed).  The  blast  wave 
direction  then  follows  along  the  direction  of  RBi  indicated  by  «,  as 


RBi  Rc+a, 
RBi  I  I  Rc+ai 


(3.15) 


The  elemental  force  acting  at  point  P  is  then  simply 


L  =  Pi(t,  |  Rc  +  «,  |  )AA(ni-ui)ui  (3.16) 


where  p,  represents  the  instantaneous  overpressure  intensity  of  the  blast  wave 
given  in  (3.1 ).  Note  that  p,  is  a  function  of  time  and  distance  between  point  P  and 
the  center  of  the  explosion. 

It  follows  that  —  assuming  the  aircraft  to  be  a  rigid  body  —  the  total  force, 
Fb,  and  moment,  MB,  exerted  by  the  blast  wave  can  be  expressed  in  terms  of  the 
sum  of  all  the  elemental  force  components  over  the  exposed  area  of  the  aircraft 
body  to  the  blast  pressure  wave: 

Fb=  Z  fi=  Z  Pi(t’\Rc+ai  DM";  •«/)«,  (3.17) 

Exposed  _  area  Exposed  _  area 

MB=  Z  aixfi=  Z  K X pt(t, I  Rc  +  \)AA(ni  ■  ui)ui ]  (3.18) 

Exposed  _  area  Exposed  _  area 


Since  momentum  must  be  conserved  in  the  process,  the  net  effect  of  the 
application  of  these  forces  over  a  time  increment  t  gives  rise  to  the  changes  of 
both  linear  and  angular  momentum  (acting  on  the  center  and  about  the  center  of 
mass)  as  follow 


F,.  -  Linear  Resistance  =  M r  — 
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(3.19) 
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dt 


J  (M B  -  Angular _Resistance)dt  =  I G  [co(t)  -  a>0  ] 


(3.20) 


where  Linear  Resistance  and  Angular  Resistance  represent  the  sum  of  all  the 
elemental  resistive  force  components.  Putting  (3.17)  and  (3.18)  into  (3.19)  and 
(3.20)  respectively  and  re-arranging  them,  we  can  express  the  variation  in  linear 
and  angular  velocity  over  time  as 
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The  presence  of  the  linear  and  angular  resistive  force  components  add  an 
additional  level  of  complexity  when  evaluating  Equations  (3.21)  and  (3.22).  In 
general,  these  components  vary  with  time  as  a  result  of  the  aircraft  motion. 
However,  even  if  the  resistive  component  effects  can  be  ignored,  the  analytical 
solution  to  (3.21)  and  (3.22)  is  only  possible  for  simple  geometrical  aircraft 
shapes  (in  that  we  also  have  to  assume  that  the  orientation  of  the  projected  area 
does  not  change  significantly  with  time  due  to  the  aircraft  motion).  Numerical 
solution  (with  the  aid  of  computer  simulation)  may  help  to  obtain  a  reasonable 
estimate  of  the  linear  and  angular  velocity  variations  for  more  general  shape  and 
different  variations.  Such  simulations  can  also  account  for  the  resistive  force 
components. 
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2.  Approach  to  Numerical  Simulation  for  ISAR  Imaging 

ISAR  imaging  of  an  aircraft  works  by  mapping  the  changes  in  the  phase  of 
the  returned  signals  to  the  range  and  cross-range  positions  of  the  target 
scattering  centers.  Hence,  our  main  interest  is  to  derive  from  (3.21)  and  (3.22) 
the  target’s  radial  range  and  angular  position  variations  that  are  related  to  the 
phases  of  the  returned  ISAR  signal  for  image  generation. 

Consider  the  following  3D  geometry  of  an  ISAR  looking  at  the  aircraft  as  it 
is  being  hit  by  the  blast  wave. 


Figure  1 1  3D  geometry  of  the  ISAR  staring  at  aircraft. 

From  the  geometry,  the  radial  velocity  detected  by  the  radar  can  be  given  by 

”'<0  =  V(')'ifoi  (3'23) 

Hence,  the  range,  R(t),  and  angular  position,  0(t)  of  the  target  over  the  entire 
ISAR  coherent  integration  period,  T,  can  be  obtained  by  integrating  (3.23)  and 
(3.22)  over  T  to  give 
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(3.24) 


tf(0  =  *o+k(0*  =  *o+f  v(0'-^7  dt 

T  T  _  I  I  _ 

0{t)  =  0O  +  J |  (o{t)  |  dt  (3.25) 

T 

From  (3.24)  and  (3.25),  we  can  obtain  the  range  variation  of  a  scattering 
point  on  the  aircraft  ( RP )  located  at  P(x,y,z)  (and  represented  by  the  vector  P )  by 
taking  the  magnitude  of  R(t)  plus  the  projection  of  the  vector  P  onto  the  fixed 
radar  frame  of  reference  (P  is  in  the  aircraft  frame  of  reference  and  its  frame  is 
rotating  with  time  varying  0  from  the  radar  frame  of  reference). 

It  is  not  trivial  to  provide  a  simple  analytical  simulation  model  for  (3.24) 
and  (3.25).  However,  if  the  distribution  of  R  and  6  can  be  obtained  over  the 
integration  period  T  and  at  intervals  equal  to  the  sampling  time  of  the  ISAR,  S„, 
through  computer-based  simulation,  the  results  of  this  simulation  can  be  used  as 
inputs  to  the  ISAR  radar  simulation  model  to  give  the  solution  for  RP  (through 
coordinate  transformation  and  projection).  Following  the  previous  discussion,  RP 
can  then  be  used  to  map  out  the  ISAR  image  through  Fourier  processing.  Figure 
12  outlines  this  methodological  approach. 


Figure  12  Approach  to  generate  of  target  motion  for  ISAR  simulation. 
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In  the  figure,  the  block  on  the  left  represents  the  simulation  model  for 
solving  the  equations  of  motion  and  provides  a  stream  of  R  and  0  data  over  the 
period  T  with  interval  to  that  of  Sn.  The  input  to  this  computer  model  includes  the 
target  parameters  (such  as  its  location  with  respect  to  the  radar,  its  shape  and 
size,  surface  orientation,  weight  and  motion),  the  blast  wave  parameters,  its 
distance  and  direction  from  the  target  center  of  mass,  the  sampling  time  and  the 
integration  period  from  the  ISAR.  The  output,  in  term  of  R  and  0  together  with  the 
target  scattering  points’  locations,  is  fed  into  the  ISAR  simulation  model  on  the 
right  where  they  are  used  to  generate  the  ISAR  image.  A  distinct  advantage  of 
using  this  approach  is  that  changes  can  be  made  to  any  one  of  the  models  with 
minimum  influence  on  the  others. 

D.  BLAST  WAVE  EFFECT  ON  AIRCRAFT  AND  ISAR  IMAGING 

In  the  development  of  Chapter  II  we  showed  that  many  ISAR  systems  rely 
heavily  on  the  linearity  of  the  target  motions  (in  radial  and  angular  velocity)  over 
the  integration  period  Tfor  radar  image  generation.  This  is  not  desirable  for  use 
of  ISAR  in  airborne  target  Battle  Damage  Assessment  (BDA)  since  the  target  is 
likely  to  possess  high  dynamic  motions  as  a  result  of  the  blast  force  interaction 
with  its  surface  body. 

In  Equation  (3.21)  and  (3.22),  we  have  shown  that  these  motions  contain 
higher  order  acceleration  components  that  will  alter  the  phases  of  the 
backscattered  radar  signal  in  a  non-linear  manner  as  it  is  collected  by  the  radar 
receiver.  Owning  to  the  finite  time  duration  of  the  blast  pressure  force,  these 
non-linear  components  only  exist  as  impulses  for  a  fraction  of  the  time  equivalent 
to  blast  duration  but  their  magnitude  can  cause  significant  changes  to  the  aircraft 
motion.  Also,  these  impulse  forces  leave  behind  residual  velocity  components 
for  the  remaining  time  of  the  radar  coherent  processing  period,  and  are  not  easily 
mitigated  by  motion  compensation  algorithms.  Hence,  both  the  presence  of  non¬ 
linear  and  residual  velocity  components  will  act  to  interfere  with  the  image 
processing  algorithm  and  distort  the  ISAR  image. 
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IV.  TIME-FREQUENCY  TRANSFORM  METHODS  OF  ISAR 

IMAGING 

A.  TIME-FREQUENCY  REPRESENTATIONS  OF  SIGNALS 

The  Fourier  transform  is  a  powerful  and  widely  used  method  in  signal 
processing  and  analysis  because  of  its  ability  to  decompose  any  arbitrary  signal 
into  a  set  of  sinusoidal  functions  with  different  frequencies.  The  distribution  of 
frequency  components  is  useful  for  characterizing  the  properties  and  behavior  of 
the  signal,  and  on  which  signal  manipulations  such  as  filtering  can  be  performed 
to  alter  them  for  specific  application  purposes.  However,  the  Fourier  operation 
and  its  inverse  involve  a  one-to-one  transitional  relationship  between  the  time 
and  frequency  domains.  Consequently,  the  “information  about  the  time 
localization  of  the  frequency  components”  cannot  be  easily  interpreted  in  the 
Fourier  spectral  domain5  (Hlawatsch  and  Boudreaux  [6]).  This  limitation  is  not  a 
problem  for  stationary  types  of  signals  but  is  generally  true  for  time-varying 
signal. 

Many  signals  encountered  in  the  real  world  have  frequency  content  that 
varies  with  time.  One  common  example  is  music  where  the  harmonic  content 
changes  for  different  notes.  Other  common  examples  include  biomedical 
signals,  speech,  vibrations  and  linear  chirp  pulses.  The  time-varying  nature  of 
such  signal’s  frequencies  makes  them  rather  difficult  to  represent  in  the  time  or 
frequency  domain  alone.  The  time-frequency  Transform  (TFT)  helps  to  alleviate 
the  situation  by  mapping  a  one-dimensional  signal  as  a  function  of  time  into  a 
two-dimensional  time-frequency  plane.  Unlike  the  Fourier  transform  (which  is 
always  linear  in  nature),  the  time-frequency  representation  of  the  signal  can  be 
linear,  bi-linear  (quadratic)  or  can  contain  higher  orders  of  non-linear  terms 
depending  on  the  method  of  signal  representation.  In  most  applications,  the  first 
two  cases  are  widely  employed  and  we  shall  discuss  them  as  follow. 


5  Strictly  speaking,  this  information  is  embedded  within  the  phase  of  the  complex  valued 
frequency  components  but  it  is  often  awkward  to  operate  on  in  the  Fourier  domain. 
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B.  SHORT  TIME-FREQUENCY  TRANSFORM  (STFT) 

1.  Definition 

The  short  time-frequency  transform  (STFT)  and  its  associated  Gabor 
transform  are  the  earliest  forms  of  linear  TFT’s  and  date  back  to  the  1940s.  In 
this  case,  the  STFT  can  be  defined  as 

/•  -jin ft 

STFT  (t,  f)  =  \s{t')w*{t'-t)e  dt'  (4.1) 

where  s(t)  is  the  signal  and  w*(t)  is  the  complex  conjugate  of  an  analysis  window 
function  with  a  finite  time  duration. 

From  (4.1),  the  STFT  can  be  seen  as  taking  the  Fourier  transform  of  the 
signal  after  it  is  multiplied  by  the  complex  conjugate  of  a  shifted  window  function 
centered  around  t.  Because  the  time  width  of  the  window  function  is  finite  and  in 
general  shorter  than  the  signal  duration,  the  effect  is  the  “suppression  of  the 
signal  outside  a  neighborhood  around  the  analysis  time  point  t’  =  t."  (Hlawatsch 
and  Boudreaux  [6])  The  STFT  represents  simply  the  local  spectrum  of  s(t’)  at  the 
analysis  time  t. 

The  magnitude  of  the  STFT(t,  f)  is  often  called  the  spectrogram  of  the 
signal  and  shows  how  the  frequency  spectrum  varies  as  a  function  of  time  in  the 
horizontal  axis.  This  magnitude  is  often  represented  by  a  surface  above  the  time- 
frequency  plane  and  displays  the  spectral  components  at  that  time. 


2.  Properties  of  STFT 

According  to  Hlawatsch  and  Boudreaux  [6],  the  STFT  also  has  a  dual 
relationship  in  the  frequency  domain  and  can  be  expressed  as 

STFT(tJ)  =  e-j2*fl\S(f')W*(f'-f)eJ2,Tr,df'  (4.2) 

where  S(f)  and  W(f)  are  the  Fourier  transform  of  s(t)  and  w(t),  respectively.  This 
relationship  gives  additional  flexibility  in  representing  the  STFT  of  the  same 
signal  from  the  spectral  domain. 


46 


There  are  other  important  properties  of  the  STFT.  The  most  significant 
property  is  that  it  satisfies  the  linearity  principle  where  the  STFT  of  the 
superposition  of  two  signals  can  be  represented  by  an  equivalent  linear 
combination  of  their  respective  signals’  STFT. 

s(t)  =  cisl(t)  +  c2s2(t) 

^STFT(t,f)  =  c1STFTl(t,f)  +  c2STFT2(t,f) 

The  linear  behavior  of  STFT  allows  the  decomposition  of  a  complex  signal  into  a 
basis  set  of  any  arbitrary  functions  to  which  the  STFT  can  be  easily  applied. 

Other  properties  include  the  preservation  of  time/frequency  under 
frequency /time  shift: 

s(t)  =  sx{t)e-j2^  =*  STFT  (t ,  /)  =  STFT^tJ  -  f0)  (4.4) 

s(t)  =  Sl(t-t0)  =>  STFTitJ )  =  STFTx{t-t0,f)e-j2xft0  (4.5) 


3.  Time  and  Frequency  Resolution 

The  time  and  frequency  resolution  of  the  STFT  are  influenced  by  the 
choice  of  window  function  and,  especially,  its  width.  From  (4.1),  good  time 
resolution  can  be  achieved  using  shorter  window  duration.  On  the  other  hand, 
the  dual  relationship  of  Equation  (4.2)  tells  us  that  a  narrow  bandwidth  is  required 
for  good  frequency  resolution.  But  the  uncertainty  principle  of  time  and 
frequency  in  the  two  domains  does  not  allow  the  existence  of  an  arbitrarily  short 
time  duration  and  small  bandwidth  window  for  the  same  signal  (Fllawatsch  and 
Boudreaux  [6]).  Hence,  this  principle  inherently  limits  the  time-frequency 
resolution  of  the  STFT. 

Consider  the  following  two  cases  (Hlawatsch  and  Boudreaux  [6]).  In  the 
first  case,  we  choose  the  window  to  be  a  Dirac  delta  function  with  perfect  time 
resolution.  Then,  for  the  signal  s(t),  its  STFT  is  given  by 

STFT(t,f )  =  J  s{t')S{t'-t)e~J2,rf,'dt'  =  s{t)ei2nf'  (4.6) 


47 


The  magnitude  is  essentially  the  same  as  the  magnitude  of  \s(t)\  without  any 
information  about  the  frequency  content  of  the  signal.  Similarly,  if  we  let  W(f)  = 
8(f)  and  apply  it  to  (4.2),  the  STFT  of  S(f)  becomes 

STFT(t,f)  =  e-12’f,\S(f)S(f-fY,Ur,df'  =  S(f)  (4.7) 

which  is  the  Fourier  transform  of  s(t)  -  its  magnitude  bears  no  information  about 
the  time  variation  of  the  frequency  components. 

The  effect  of  time  and  frequency  resolution  can  also  be  seen  by  referring 
to  the  illustrative  example  as  follows.  Figure  14(a)  shows  a  signal  whose 
frequency  changes  with  time  and  Figure  14(b)  is  the  Fourier  transform  of  this 
signal.  From  Figure  14(b),  we  can  see  that  the  frequency  domain  does  not 
provide  any  indication  on  how  the  signal  frequency  changes  with  time. 


(a)  Time  Domain  Signal 
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(b)  Frequency  Domain  Representation 


Figure  13  Example  of  a  time-varying  signal  with  changing  frequency. 
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The  time-frequency  representation  of  the  same  signal  is  shown  in  Figure 
14.  This  figure  was  generated  using  a  sliding  Gaussian  window  function.  In 
Case  1 ,  we  set  the  width  of  the  window  function  such  that  good  time  resolution 
can  be  obtained:  notice  the  significant  overlap  in  frequency.  In  Case  2,  we 
broaden  the  window  width  to  give  good  frequency  resolution  but  in  turn  the  time 
resolution  suffers. 
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(a)  Case  1 :  Small  window  size  for  good  time  resolution. 
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(b)  Case  2:  Large  window  size  for  good  frequency  resolution. 


Figure  14  Effect  of  Time  and  frequency  resolution. 
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4.  Shape  of  the  Window  Function 

The  specific  shape  of  the  window  function  is  important  in  reducing  the 
side-lobe  interference  in  the  STFT  representation.  Typically,  one  requires  the 
window  function  to  taper  toward  zero  smoothly  and  applies  the  Hamming, 
Hanning,  Kiaser-Bessel,  or  Gaussian  window  function  (Chen  [2]).  The  best 
results  have  been  shown  to  be  achieved  by  Gaussian  windows  where  the  time- 
bandwidth  resolution  product  is  AtAf=l/2  (Refer  to  Chen  [2]).  The  use  of 
Gaussian  windows  in  an  STFT  is  known  as  a  Gabor  transform. 


C.  BI-LINEAR  TIME-FREQUENCY  TRANSFORM 

Although  the  STFT  form  of  signal  representations  has  a  nice  linear 
property,  the  time-frequency  resolution  is  limited  by  the  uncertainty  principle, 
which  at  best  gives  a  time-bandwidth  product  of  1/2.  Accordingly,  Cohen  [12], 
Hlawatsch  and  Boudreaux  [6]  and  many  others  have  shown  that  significant 
improvement  in  the  time-frequency  resolution  can  be  obtained  by  representing 
the  signal  in  the  quadratic  or  bi-linear  TFT  form.  The  quadratic  representation  is 
also  “an  intuitively  reasonable  assumption  when  we  want  to  interpret”  the  time- 
frequency  of  the  signal  as  a  form  of  energy  distribution  (Hlawatsch  and 
Boudreaux  [6]). 

To  classify  a  TFT  signal  representation  in  the  bi-linear  form,  it  must  satisfy 
the  following  marginal  properties  related  to  the  energy  density. 


\TFT(t,f)df  =  p{t)  =  \s(t)\2 
\TFT(t,f)dt  =  P(f)  =  \S{f)\2 


(4.8) 


where  p(t)  and  P(f)  are  the  instantaneous  power  and  the  spectral  energy  density 
of  the  signal  s(t)  respectively.  As  a  consequence,  the  signal  energy  E  =  I  \x(t)\2dt 
=  I  \X(t)\2df  can  be  derived  by  integrating  TFT(t,f)  over  the  entire  time-frequency 
plane. 
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The  most  difficult  part  in  interpreting  the  bi-linear  TFT  representations  is 
the  presence  of  the  cross-interference  terms  which  occur  as  a  result  of  the 
quadratic  superposition  principle.  Consider  the  spectrogram  of  the  STFT  (which 
is  often  loosely  interpreted  as  a  bi-linear  TFT  representation  although  it  does  not 
meet  the  marginal  properties  according  to  Hlawatsch  and  Boudreaux  [6]). 

SPEC(t,f)  =  \STFT(t,f)\2  (4.9) 

We  can  easily  show  that,  for  the  signal  s(t)  =  si(t)  +  s2(t),  the  spectrogram  of  this 
composite  signal  is  not  SPECi(t,f)  +  SPEC2(t,f)  but  contains  cross  interference 
terms: 


SPEC(t,  /)  =  |  STFT,  (t,  /)  +  STFT2  (t,  f)f 

=  \STFTt  (t,  /)  |2  +  \STFT2  ( t ,  /)  |2  +  2STFT,  (t,  f)STFT2  (t,  /)  (4.10) 

=  SPECl  (t,  f )  +  SPEC 2  (/,/)  +  Cross  _  Intereference 


This  behavior  is  generally  true  for  all  types  of  bilinear  TFT  representations 
although  the  magnitude  of  the  cross  term  varies  with  different  types  of  bi-linear 
TFT  transforms.  Also,  the  more  complex  the  signal  is,  the  more  the  linear 
combination  of  its  basis  set  is  used  to  represent  the  signal;  hence,  the  more 
cross  interference  terms  that  will  appear  in  the  spectrogram  and  bi-linear  TFT 
representations  (Hlawatsch  and  Boudreaux  [6]). 


There  are  many  forms  of  bilinear  TFT.  The  most  basic  is  the  Wigner-Ville 
Distribution  (WVD),  which  is  defined  as  the  Fourier  transform  of  the  time- 
dependent  autocorrelation  of  the  signal  R(t,  t)  (Chen  [2]) 

WVD(t,  /)  =  J  R(t,t  ')e-j2*f,'dt  ’  (4.11) 


where  R(t,  t)  is  chosen  as 


R(t,t')  =  s 


f  f) 

s 

\  t  H - 

s  * 

l  2  J 

l  2  J 

(4.12) 


Because  of  this  relationship  with  the  auto-correlation  function,  the  WVD  is 
known  to  have  very  good  time-frequency  localization  and  hence  good  time  and 
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frequency  resolution.  But  this  distribution  is  also  known  to  possess  large  cross 
interference  terms  that  hinder  its  direct  usefulness  for  signal  analysis  and 
processing. 

A  more  generalized  form  known  as  the  “Cohen  class”  of  transforms 
extends  the  use  of  the  WVD  to  include  other  members  of  the  bilinear  TFT 
representation  with  a  kernel  function  <j)(t,  t).  The  general  form  of  the  Cohen  class 
is  given  by  (Chen  [2]) 


f  t" 
l  2y 


V 

u - 

.  2  j 


<j>(t-u,t')e  ilnf>  dudt  ’ 


(4.13) 


Note  that  the  WVD  is  simply  a  subset  of  the  Cohen’s  class  with  <f>(t,  t)  =  S(t).  The 
significance  of  the  Cohen  class  is  that  different  types  of  kernel  function  can  be 
designed  to  reduce  the  cross-term  interference  of  the  bilinear  TFT  representation 
but,  in  the  process,  compromise  the  time-frequency  resolution  (compared  to 
WVD  —  which  is  known  (Chen  [2])  to  possess  the  best  time-frequency 
resolution).  Two  well-known  distributions  in  this  category  are  the  Choi-Williams 
distribution  (CWD)  and  the  cone-shaped  distribution  (CSD).  Further 
development  of  these  topics  is  beyond  the  scope  of  this  thesis,  and  we  shall  not 
discuss  them  further.  For  detailed  reviews,  refer  to  excellent  texts  by  Chen  [2], 
Fllawatsch  and  Boudreaux  [6]  and  Cohen  [12], 


D.  TIME-FREQUENCY  BASED  IMAGE  FORMATION 

In  Chapter  II,  we  noted  the  stringent  radial  and  angular  motion 
requirements  imposed  on  ISAR  systems  for  high  quality  image  generation:  both 
velocities  must  essentially  remain  constant  throughout  the  coherent  integration 
period.  The  constant  radial  velocity  allows  the  radar  to  accurately  estimate  the 
target’s  linear  motion  and  compensate  for  it  accordingly.  To  maintain  a  simple 
direct  relationship  for  mapping  the  returned  phase  variations  to  the  cross-range 
position  of  the  target’s  scattering  centers,  the  angular  rotation  rate  must  be  linear 
and  within  an  upper  bound  such  that  the  viewing  angle  of  the  ISAR  is  small 
enough  to  prevent  cell  migration. 
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However,  when  the  aircraft  is  being  hit  by  explosive  munitions,  the  blast 
force  will  cause  changes  in  both  linear  and  angular  momentums.  The  net  effect 
is  that  the  target  will  exhibit  complex  variations  in  its  radial  and  angular  velocities, 
occurring  over  a  very  short  time  interval  (typically,  msec)  and  well  within  the 
coherent  integration  time  (usually  measured  in  sec).  Most  translational  and 
rotational  motion  compensation  techniques  are  inadequate  to  provide  good  range 
and  Doppler  tracking  of  these  complex  motions.  As  a  result,  the  image  becomes 
distorted  when  conventional  Fourier-based  techniques  are  employed  for  image 
construction. 

The  complex  variation  of  the  radial  and  angular  velocities  from  the  blast  hit 
can  be  viewed  as  a  time-varying  frequency  shift  of  the  returned  signals  within  the 
integration  period.  Although  these  signals  pose  problems  for  conventional  FFT 
processing  in  ISAR  systems,  the  signals  do  contain  information  about  the 
dynamic  behavior  of  the  aircraft  as  it  is  hit,  and  this  information  can  be  extracted 
using  TFT  processing. 

From  an  ISAR  system  architecture  point  of  view,  the  front  end  processing 
remains  largely  the  same.  The  SFW  returned  signals  are  detected  and  sampled 
as  I-  and  Q-  data,  and  shipped  off  to  the  appropriate  range  and  Doppler  tracking 
algorithm  for  motion  compensation.  Synthetic  range  profile  processing  is  then 
performed  on  the  M x  N  complex-valued  sampled  data  array  to  obtain  the  range 
information  as  usual.  The  result  is  also  a  M xN  complex-valued  matrix  array,  H. 
Note  that,  in  accordance  with  our  earlier  definition,  each  row  of  the  matrix 
(denoted  as  h(  m  ,  1:N )  where  m  =  1,  2,  M  )  represents  a  time-history  series 
comprised  of  N  pulses  of  range  cells.  Each  column  of  the  matrix  (denoted  as  h( 
1:M ,  n  )  where  n  =  1,  2,  N)  represents  a  collection  of  time-histories  of  the 
target  from  M  burst  at  the  particular  range  cell  index  n. 

In  conventional  ISAR  image  processing,  the  direct  Fourier  transform 
(usually  implemented  by  an  FFT)  is  applied  to  each  column  h(  1:M ,  n  )  to  obtain 
the  cross-range  position  of  the  target  scattering  points  in  that  range  cell  index. 
The  Fourier  processing  still  yields  a  one  dimensional  column  vector: 
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DFFT(\:M  ,n)  =  F FT n  {h(l:M,n)}  (4.14) 

This  process  is  repeated  for  all  N  range  cell  indices.  The  result  is  still  a  2D  MxN 
matrix  array  with  the  magnitude  of  each  element  indicating  the  presence  of  the 
target  point  scatterer  in  the  particular  range  and  cross-range  positions. 

TFT’s  differ  from  the  conventional  Fourier  transform  method  of  processing 
in  that  the  windowed  transform  (can  be  linear  STFT  or  bilinear  TFT)  maps  the 
one-dimensional  collection  of  time-histories  of  M  bursts  in  each  column  h( : ,  n  ) 
into  a  two-dimensional  M  x  M  plane  revealing  information  about  the  Doppler 
variation  along  the  time-histories  in  that  cell  index. 

Dtft (1 :  M,  1 : M, n)  =  TFTn  {h{  1 : M ,nj)  (4.15) 

This  process  is  repeated  for  all  range  cells  and  we  get  a  three  dimensional 
complex-valued  array  of  Mx  Mx  N  Doppler-time-range  matrix.  Figure  15  depicts 
the  image  processing  based  on  TFT  method. 
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Figure  15  ISAR  image  processing  based  of  TFT  method. 


A  series  of  M  ISAR  images  can  then  be  obtained  from  the  three 
dimensional  matrix  array  by  summing  over  all  the  range  cell  indices  in  each  m 
time-history  point  as  shown  in  the  figure  below. 
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M  Doppler 
Mapping 
from  TFT 


Figure  16  Extraction  of  ISAR  Image  from  the  Doppler-time-Range  Matrix. 


The  image  extraction  process  can  be  summarized  by  the  following  expression: 

Im (1 :  M, l:N)  =  D(\ :  M, m,  1 :  N)  for  m  =  1,  2,  ....  M  (4.16) 

When  the  same  process  is  repeated  for  each  time-history  series,  the  time-history 
series  display  of  the  2D  ISAR  image  with  indices  running  from  m  =  1  to  M  is 
analogous  to  a  moving  picture  of  the  aircraft  within  the  integration  period  (as 
opposed  to  the  static  one  by  the  conventional  Fourier-based  construction 
method).  Hence  the  dynamic  motion  of  the  aircraft  due  to  the  blast  effect  is 
revealed  by  these  sets  of  ISAR  images. 

The  image  resolution  from  TFT’s  is  generally  coarser  in  comparison  to  the 
image  obtained  by  conventional  FFT.  The  resolution  also  depends  on  the  width 
and  shape  of  the  window  function  for  linear  TFT.  Bilinear  TFT’s  will  give  better 
resolution,  but  the  presence  of  cross  interference  terms  may  act  to  distort  the 
image. 
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V.  DEMONSTRATION  OF  BLAST  EFFECT 
CAPTURINGBY  ISAR  USING  TIME-FREQUENCY 

TRANSFORM 

A.  BLAST  EFFECTS  AND  THE  ISAR  SIMULATION  MODEL 

To  demonstrate  the  possibility  of  using  Time-Frequency  Transform  (TFT) 
methods  to  ‘see’  the  dynamic  behavior  of  an  aircraft  target  in  ISAR  imaging,  a 
simple  simulation  model  was  developed  as  part  of  this  thesis.  This  simulation 
contains  two  primary  sections  following  the  methodological  approach  described 
in  Chapter  III.  The  first  section  aims  to  simulate  the  response  of  the  aircraft  as 
the  result  of  a  blast  hit.  The  second  part  of  the  model  extends  the  usual  ISAR 
imaging  simulation  model6  to  include  a  time-frequency  method  of  radar  imaging 
based  on  the  work  of  Chen  [2], 

The  simple  model  developed  for  this  demonstration  does  not  provide  a 
complete  description  of  the  real  world  problem  at  hand  involving  target  BDA. 
There  are  several  limitations  associated  with  this  model  that  deserve  comment: 

1)  Only  two-dimensional  geometries  are  considered  to  ease  the 
development  and  thesis  illustration.  Nevertheless,  one  can  extend  the 
model  to  the  three-dimensional  case  with  careful  considerations  about  the 
target  roll,  pitch  and  yaw  effects,  and  with  some  manipulations  of  the 
frame  of  references  using  matrix  transformation  and  vector  operations. 

2)  Initially,  we  shall  consider  only  the  translational  response  of  the  air 
target  in  this  demonstration.  The  model’s  implementation  requires 
analytical  expressions  for  translational  and  rotational  velocities.  It  is 
easier  to  obtain  the  translational  motions  from  (3.21)  with  some 
reasonable  assumptions  (as  we  shall  see  later).  However,  the  angular 
dependence  in  Equation  (3.22)  contains  terms  related  to  the  geometrical 

6  Some  of  these  models  are  easily  available  in  or  as  part  of  publications  related  to  ISAR 
imaging  in  the  courtesy  of  the  author  (s).  For  this  thesis,  we  would  like  to  acknowledge  that  our 
model  is  a  modified  version  of  the  Matlab  ISAR  simulation  program  using  conventional  FFT 
processing  from  Jae,  Thomas  and  Flores  [3]  to  help  in  hastening  the  model  development  time 
compared  to  starting  from  scratch. 
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distribution  of  the  projected  aircraft  area  and  the  distribution  of  the 
pressure  wave  along  its  surface  that  needs  to  be  properly  accounted  for  in 
order  to  obtain  reasonable  estimates  of  the  net  effect  in  angular  velocities. 
There  is  no  simple  analytical  solution  to  this  problem,  and  it  requires  a 
large  scale  numerical  approach  which  is  beyond  the  scope  of  this  thesis. 

3)  Linear  time-frequency  transform  methods  (in  particular  the  Short 
Time-Frequency  Transform  or  STFT)  will  be  used  in  this  study.  The  main 
concern  with  linear  TFT’s  is  the  time-frequency  resolution,  which  (at  best) 
gives  a  time-bandwidth  product  of  only  !4  We  will  have  to  evaluate  its 
adequacy  for  our  purpose  before  we  proposed  further  improvements  to  it. 

The  code  for  the  simulation  program  was  developed  using  Matlab  6.5  with 
time-frequency  toolbox  provided  for  by  Auger  et  al  [7]  &  [8]  available  at 
http://crttsn.univ-nantes.fr/~auqer/tftb.html,  Aug  2004.  The  entire  code  is  given  in 
Appendix  A. 

B.  SIMULATION  OF  AIRCRAFT  RESPONSE  TO  BLAST  WAVE 
INTERACTION 

The  primary  goal  of  this  part  of  the  simulation  is  to  obtain  the  translational 
motion  response  of  the  aircraft  from  the  blast  hit.  This  response  appears  as  a 
form  of  variation  in  range  as  a  function  of  time  and  can  be  derived  from  the 
equation  of  motion  for  translational  motion  (3.21)  and  the  blast  wave  expressions 
given  in  (3.1).  There  are  certain  assumptions  needed  in  order  to  simplify  the 
derivation.  They  are  described  as  follows. 

1.  Blast  Pressure  Distribution  on  Aircraft 

To  ease  the  stringent  geometrical  requirement  of  the  non-symmetrical 
aircraft  surface,  we  assume  the  blast  wave  to  be  a  plane  wave  with  equal 
pressure  distribution  on  the  exposed  area  of  the  aircraft  body.  The  exposed  area 
is  also  assumed  to  remain  constant  with  time  during  the  evaluation  period.  By 
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ignoring  the  resistive  force  components,  Equation  (3.21)  can  be  re-arranged  and 
simplified  to  become 


v(t)  =  v„+-^-[\p(t>Rc)dt  Z  AA 


M, 
A , 


G  \t 


(«,  •  «,)«,• 


=  V  +- 


M, 
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J  p(t,Rc)dt 

K  t  J 


J  y  Exposed  _  area  J 
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(5.1) 


where  v„  is  the  initial  velocity,  MG  is  the  mass  of  the  aircraft,  Ap  is  the  projected 
area  of  the  aircraft  exposed  to  the  blast  pressure  given  by 


A,= 


X  M 

Exposed  _  area  J 


(«/•«;) 


(5.2) 


and  p(t,  Rc)  is  the  instantaneous  pressure  wave  described  by  Equation  (3.1 ).  Rc  is 
factored  into  the  original  p(t)  since  the  parameters  determining  the  shape  of  p(t) 
—  Po,  ta,  td  and  a  —  are  actually  dependent  on  the  blast  distance,  Rc.  The  unit 
vector  Hi  represents  the  direction  of  the  net  blast  effect  which  is  normal  to  the 
aircraft  surface  subjected  to  the  blast. 


Now  that  the  blast  pressure  acting  on  the  aircraft  body  is  independent  of 
its  direction  and  that  the  distance  Rc  can  be  taken  to  have  constant  value  under 
the  planar  wave  assumption,  the  impulse  response  (which  is  integral  of  p(t)  over 
time  t )  can  be  evaluated  separately.  If  we  denote  the  impulse  response  of  the 
target  by  I(t),  then 


I(t)  =  \p(t',Rc)dt' 


1- 


t'-taW 

td(K) 


exp 


td{Rc )  , 


(5.3) 


dt '  where  t<t<  t,. 


where  p0(Rc),  ta(Rc),  td(Rc)  and  a(Rc)  can  be  expressed  as  a  function  of  the  blast 
distance  using  the  empirical  Equations  (3.2),  (3.4),  (3.6)  and  (3.8)  respectively. 

We  denote pa’  as p0(Rc),  ta’  as  ta(Rc),  td’  as  td(Rc)  and  a’  as  a(Rc).  Applying 
the  appropriate  boundary  conditions,  we  can  obtain  the  solution  fo r  I(t)  as 
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2.  Geometry  of  the  Problem  and  its  Analytical  Solution 

We  are  interested  in  the  radial  velocity  component  since  this  is  detectable 
by  the  radar  as  a  Doppler  shift.  To  obtain  the  radial  velocity  variation,  consider 
an  exaggerated  two-dimensional  geometry  of  an  ISAR  looking  at  the  aircraft 
while  an  explosive  blast  occurs  near  it  as  shown  in  Figure  17. 


Figure  17  2D  Geometry  between  ISAR,  aircraft  and  blast  point. 

For  the  geometry  of  the  figure,  the  force  acting  along  the  radar  line-of- 
sight  (LOS)  is  given  by 

FLOS=ApP(1’RC)Sin0B  (5-5) 
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where  Ap  is  the  projected  surface  area  of  the  aircraft  to  the  blast,  p(t,  Rc)  is  the 
instantaneous  blast  pressure  wave  as  in  (3.1),  Rc  is  the  blast  distance,  and  0B  is 
the  direction  of  the  blast  wave. 

We  use  the  following  sign  convention:  the  radial  velocity  is  positive  when 
the  target  is  approaching  the  radar  and  negative  otherwise.  Hence,  the  resultant 
radial  velocity  variation  can  be  simply  obtained  using  conservation  of  momentum 
(ignoring  all  other  resistive  forces)  as  follows 

-J Fws (*)&  =  Mg  [vr(f)- vr(0)]  (5.6) 


Re-arranging  (5.6)  and  substituting  (5.5)  into  it,  we  get 


An  sin  6 R  r 

K(0  =  v,-(°) - ~TZ - J  p(t,Rc)dt 


M, 


G  I 


An  sin  6 R 

=  Vr(0) - — - ~ 
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m 


(Equation  (5.7)  can  also  be  obtained  from  the  projection  of  (5.1)  onto  the  radar 
LOS  and  the  integral  is  exactly  the  same  as  in  Equation  (5.4).) 

For  the  ISAR  simulation  model,  the  target  linear  motion  is  expressed  as 
the  range  variation  R(t).  To  obtain  R(t),  we  simply  integrate  vr(t)  with  respect  to 
time  using  the  following  kinematic  relationships. 

vXt')=Cy-^\dR  =  \  vr  (t  ')dt ' 
at  R  t 

R(t)  =  R(Q)  +  \vXt')dt'  (5.8) 

t 

An  sin  0R  r 

=  *(0)  +  v,(0)f-  p  B\l{t')dt' 

t 


The  result  of  the  integral  I(t)  is  expressed  below,  which  represents  the 
translational  dynamic  behavior  of  the  aircraft  due  to  blast  effect. 
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3.  Aircraft  Profile 

Following  Equation  (5.2),  we  can  define  the  overall  projected  area  of  the 
aircraft  exposed  to  the  blast  as  the  summation  of  the  projected  frontal  area  of  its 
body,  Af,  and  the  projected  side  area  of  its  body,  As,  and  expressed  in  term  of  the 
aircraft’s  orientation  with  respect  to  the  blast  direction  as  in  Figure  17. 

Ap  =  | Af  cos {6b  +  6) |  + 1 As  sin(6^  +  6) |  (5.11) 

The  front  and  side  profiles  of  the  aircraft  used  in  this  simulation  have  a 
shape  shown  in  Figure  18.  Hence,  AF  and  As  are  given  by 

h2 

Af  ~  n  +  (^5  —  h)hw  (5  12) 

As  =  hL  +  LThT 
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Figure  18  Aircraft  area  profile. 

4.  Blast  Parameter 

Peak  overpressure  p„,  arrival  time  ta,  blast  duration  td  and  the  waveform 
parameter  a  are  needed  to  fully  describe  the  blast  wave  as  a  function  of  time.  pQ, 
td  and  a  can  be  easily  obtained  using  the  empirical  equations  given  in  Chapter  III 
—  (3.2),  (3.6)  and  (3.8)  respectively.  The  arrival  time  ta,  however,  requires 
integration  of  the  speed  from  the  center  of  explosion  to  the  distance  the  shock 
front  hits  the  aircraft.  To  simplify  the  simulation,  tabulated  forms  of  these  data 
are  extracted  from  Kinney  and  Graham  [10]  where  they  are  scaled  to  lkg  TNT 
equivalent  values  at  standard  atmospheric  pressure  and  temperature  conditions. 
The  actual  values  can  be  obtained  by  using  Hopkinson’s  scaling  law  described  in 
Chapter  III. 

For  our  simulation,  the  values  for  the  explosive  parameters  are:  explosive 
weight  =  35kg  with  a  TNT  equivalent  ratio  of  1 .5  for  HE  type  of  explosive.  These 
are  figures  appropriate  to  many  long  range  Surface-to-air  and  Air-to-air  missiles. 
The  blast  distance  is  chosen  to  be  5m,  which  is  within  the  miss  distance  of  most 
missiles  today. 
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C.  SIMULATION  OF  ISAR  IMAGING  USING  TIME-FREQUENCY 

TRANSFORM 

The  simulation  model  for  generation  of  the  ISAR  image  consists  of  three 
distinct  parts:  (1)  Target  model  representation,  (2)  Radar  signal  generation  and 
(3)  ISAR  Image  processor.  The  description  of  each  of  these  portions  is  as 
follows. 

1.  Target  Model  Representation 

It  is  a  common  practice  in  most  ISAR  simulations  to  represent  the  ISAR 
target  as  a  series  of  point  scatterers  outlining  the  shape  of  the  target.  The 
number  of  scattering  points  and  the  reflectivity  of  each  point  can  be  arbitrary 
chosen  as  long  as  they  are  within  the  image  plane.  For  our  case,  49  scattering 
points  were  used  to  outline  a  target  and  their  reflectivity  set  to  one  (i.e.,  uniform 
signal  reflection). 

In  the  model,  the  target  is  free  to  assume  any  initial  orientation  with 
respect  to  the  radar  by  varying  60.  The  motion  for  each  of  the  target  scattering 
points  is  modeled  in  term  of  range  variation  using  Equation  (2.4),  and  since  the 
target  is  treated  as  a  rigid  body,  the  range  variation  for  all  the  points  is  influenced 
by  the  same  kinematic  equations  for  R(t)  and  G(t).  The  expression  used  for  R(t) 
follows  from  the  discussion  of  section  (5.8),  and  6(t)  is  assumed  to  behave  as  in 
Equation  (2.3). 

2.  Radar  Signal  Generation 

The  image  plane  for  the  ISAR  processing  is  represented  by  64  x  64  cell 
indices  —  a  N  x  M  matrix  array  (note  that  the  range  and  cross-range  axis  are 
swapped  in  this  presentation).  The  radar  parameters,  fc,  Af,  N  pulses,  M  burst, 
pulse  interval  1/  PRF  and  integration  period  T  together  with  the  initial  range  R0 
and  rotation  rate  co0  are  specifically  chosen  to  optimize  the  overall  radar  system 
such  that  an  ideal  image  can  be  obtained  without  the  need  to  consider  range 
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alignment,  range  and  Doppler  tracking,  and  motion  compensation  in  the  process. 
A  summary  of  these  parameters  are  given  in  Table  1  below. 

The  radar  signal  generation  follows  the  theory  in  Chapter  II  and  is  simply: 
compute  the  radial  range  R(t),  angular  variation  0(t)  and  stepped  frequency 
increases  fnM  under  (2.8)  at  each  sampling  time  tm,n  in  (2.9)  for  T,  and  re-arrange 
them  as  a  N  x  M  matrix.  These  three  time-varying  matrices  are  imposed  as 
phases  onto  the  target  model  and  the  signal  is  taken  as  the  sum  of  these  phases 
from  all  the  scattering  points  according  to  (2.10). 


Table  1 .  Radar  Parameters 


Parameters 

Values 

Remarks 

Target  Initial  Rotation  Rate,  co 

3.77deg/s 

Optimized  for  the  target  length. 

Target  Initial  Range,  Ro 

22.567km 

Optimized  for  range  alignment. 

Carrier  Frequency,^ 

3GHz 

Typical  Operating  value. 

Stepped  Frequency,  Af 

2.2MHz 

Optimized  for  range  resolution. 

Number  of  Pulses,  N 

64 

Arbitrary  chosen. 

Number  of  Burst,  M 

64 

Arbitrary  chosen. 

Pulse  Interval,  T2 

150psec 

Optimized  for  Doppler  resolution. 

Pulse  Width,  T[ 

0.454psec 

Tj=l/Af 

Integration  Period,  T 

0.6144sec 

T  =  N xMx  T2 

Max  Target  Length,  x 

30m 

x=2(Arc)~/A 

Range  Resolution,  Ar 

1.06m 

Ar=c/(2NAf) 

Cross-Range  Resolution,  Arc 

1.23m 

Arc=A/(2coT) 
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3.  ISAR  Image  Processor 

The  image  processing  of  the  radar  signal  begins  by  generating  the 
synthetic  range  profile  of  the  target  according  to  Equation  (2.1 1).  Subsequently, 
the  resultant  data  is  passed  to  another  FFT  process  to  construct  the  conventional 
two-dimensional  mapping  of  the  target  image  using  (2.12). 

On  the  other  hand,  time-frequency  image  construction  follows  the 
procedures  discussed  in  Chapter  IV:  the  STFT  is  applied  to  each  row  of  the  Nx 
M  synthetic  range  profile  matrix  as  in  Equation  (4.15).  The  result  is  a  three- 
dimensional  Doppler-time-range  mapping  of  the  target.  Each  of  the  M  time- 
histories  is  then  extracted  for  display  accordingly  as  shown  in  Figure  16.  The 
operation  is  carried  out  using  Equation  (4.16). 
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VI.  RESULTS  AND  DISCUSSIONS 


A.  TIME-FREQUENCY  REPRESENTATION  OF  ISAR  IMAGE 
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(a)  Target  model.  (b)  ISAR  image. 

Figure  19  Target  model  representation  and  its  ISAR  Image  (from  FFT  process). 

Figure  19(a)  above  show  the  model  of  an  aircraft  used  in  this  study.  The 
model  is  represented  by  the  series  of  49  point  radar  scatterers  spread  over  the 
body  length  of  35m  and  wing  span  of  20m.  Figure  19(b)  is  an  ISAR  image  of  the 
aircraft  model  generated  using  the  conventional  Fourier  process  (in  particular 
Fast  Fourier  Transform,  FFT)  described  in  Chapter  II. 

When  the  same  target  model  is  analyzed  using  the  Short  Time-Frequency 
Transform  (STFT)  method  (discussed  in  Chapter  IV),  a  successive  series  of  64 
time-history  ISAR  images  can  be  displayed  instead.  Each  of  these  images  shows 
the  same  outline  of  the  aircraft  image.  Figure  20  below  shows  an  example  of  16 
images  extracted  from  the  time-history  series  of  m  =  33  to  48  (at  time  t  between 
0.3168  to  0.4608sec  where  t  =  0  indicates  the  beginning  of  the  coherent  processing 
period). 

The  image  resolution  is  comparatively  poorer  —  because  of  the  tradeoff  in 

frequency  resolution  given  to  the  time  resolution.  Conventional  FFT’s  operate 

over  the  entire  time-history  series  of  M  cells  (in  our  case  M  =  64)  but  produce 
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only  a  one-dimensional  mapping  in  the  spectral  domain.  Whereas  in  STFT,  the 
Fourier  processing  is  carried  out  with  a  sliding  time  window  function  to  yield  a 
two-dimensional  time-frequency  mapping  of  each  time-history  series  m  where  m 
=  1,  2,  M.  The  width  of  the  window  is  generally  chosen  to  be  shorter  than  Mto 
give  better  time-resolution,  which  is  important  (as  we  shall  see  later)  for 
observing  the  dynamic  behavior  of  the  target  within  the  processing  period.  (In  its 
extreme,  we  can  show  that  exact  FFT  image  is  obtainable  from  the  STFT  simply 
by  setting  the  window  width  to  cover  the  entire  M  cells).  But  the  instantaneous 
frequency  resolution  suffers  as  a  result,  and  this  frequency  is  directly  related  to 
the  cross-range  resolution  of  the  ISAR  image. 
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Figure  20  Sample  of  ISAR  from  STFT. 
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B.  TRANSLATIONAL  RESPONSE  OF  AIRCRAFT  DUE  TO  BLAST 

EFFECT 

1.  Blast  Characteristic 

In  this  study,  we  place  the  blast  directly  between  the  target  and  radar 
along  its  line-of-sight  (LOS)  so  that  the  magnitude  of  translational  response  will 
be  maximal  in  the  radial  direction.  In  the  first  case  studied,  we  only  let  the  side 
of  the  aircraft  be  exposed  to  the  blast  wave.  The  aircraft  has  a  side  area  of  64m2 
and  weight  of  22tons.  The  translational  responses  in  term  of  acceleration, 
velocity  and  range  variation  are  as  shown  in  Figure  21  and  Figure  22. 


Figure  21  Target  acceleration  and  velocity  variation  (Blown  up  version). 


Figure  22  Target  range  variation. 
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As  expected,  the  overall  acceleration  directly  follows  the  impulse  response 
of  the  blast  pressure  wave.  The  presence  of  this  acceleration  force  gives  rise  to 
the  non-linear  variation  in  the  radial  velocity  which  eventually  dissipates  away.  By 
conservation  of  momentum,  the  change  in  velocity  will  cease  as  the  force  is 
taken  off  the  load  and  is  set  to  a  new  (constant)  value.  We  note  also  that  the 
range  variation  from  zero  to  3m  results  from  the  velocity  components.  Since 
these  effects  happen  within  one  integration  period,  the  ISAR  system  will  not  have 
enough  time  to  compensate  for  them.  Hence  the  blast  effects  are  effectively 
introduced  as  errors  in  the  imaging  process. 

2.  ISAR  Image  from  Conventional  FFT 

These  errors  are  manifest  as  distortions  to  the  radar  image  as  shown  in 
Figure  23(b)  when  constructed  using  a  conventional  FFT.  Notice  the  ‘smearing’ 
of  the  scattering  points  due  to  the  migration  of  the  range  and  cross-range  cells  in 
the  presence  of  the  radial  velocity  component,  and  the  creation  of  a  seemingly 
‘double  image.’  We  shall  discuss  this  later  when  we  examine  the  STFT  image. 
But  other  than  the  presence  of  cell  migration,  the  ISAR  image  does  not  reveal 
anything  about  the  target’s  behavior  within  the  integration  time. 


10  20  30  40  50  60 

Cross-Range  Bin 


(a)  Undistorted.  (b)  Distortion  due  to  blast  effect. 

Figure  23  ISAR  image. 
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3.  ISAR  Image  from  STFT 

The  ISAR  images  constructed  using  STFT  are  shown  below  in  Figure  24. 
The  period  of  interest  is  around  the  time  when  the  blast  interaction  occurs,  which 
is  the  first  0.3sec.  In  the  images,  we  can  observe  the  sequence  of  events  that 
happen  to  the  aircraft  when  it  is  being  hit  by  the  blast  wave.  During  the  first 
0.15sec,  the  position  of  aircraft  image  is  the  same  as  before.  As  the  blast  wave 
arrives,  the  aircraft  image  begins  to  blur  because  of  the  velocity  variations.  This 
happens  for  about  the  same  duration  as  the  blast  wave  interaction  with  the 
aircraft  body.  When  the  variation  stops,  the  aircraft  image  becomes  gradually 
clearer  (at  t  =  0.2sec)  but  seems  to  be  offset  to  a  new  position,  which  is  an 
expected  effect  of  a  result  of  the  uncompensated  velocity  components  present  in 
the  image  processing. 


TH1  (t  =  0  0096sec) 


TH5  (t  =  0  048sec) 


TH2  (t  =  0.0192sec) 


TH6  (t  =  0  0576sec) 


TH10  (t  =  0.096sec) 


TH14  (t  =  0.1344ssc) 


TH3  (t  =  0  0288sec) 


TH7  (t  =  0  0672sec) 


TH11  (t  =  0.1056sec) 


TH15  (t  =  0-144sec) 


TH4  (t  =  0  0384sec) 


TH8  (t  =  0  0768sec) 


TH12  (t  =  0.1152sec) 


TH16  (t  =  0  1536sec) 


(a)  At  time  t  =  0  to  0.1536sec. 
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(b)  At  time  t  =  0.1536  to  0.3072sec. 


Figure  24  STFT  ISAR  images  for  observing  the  blast  effect. 

From  these  STFT  images,  we  also  make  two  interesting  observations: 

1 )  There  is  a  direct  correlation  with  these  images  to  the  one  obtained 
using  FFT.  We  notice  that  the  ‘smearing’  of  the  image  earlier  is  actually 
caused  by  the  ‘excessive’  cell  migration  from  the  nonlinear  variation  of  the 
radial  velocity  when  the  blast  wave  hits  the  aircraft.  Also,  the  brighter 
‘double  image’  is  simply  the  result  of  continuous  imposing  of  the  offset 
aircraft  image  after  the  blast.  However,  note  that  the  blast  has  unusually 
large  effect  along  the  cross-range  axis  rather  than  the  range  (The  cross¬ 
range  cell  shift  is  counted  to  be  about  8-10  cells  whereas  the  range  cell 
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variation  is  only  3  cells)  although  the  only  velocity  component  is  presence 
in  the  radial  direction. 

2)  The  image  resolution  is  much  poorer  compared  to  Figure  20.  This 
is  a  consequence  of  reducing  the  size  of  the  window  function  to  only  !4  of 
the  M  cells  so  that  the  time  resolution  is  appropriate  to  reveal  the  blast 
event  happening  around  the  time  t  =  0.165  to  0.19sec.  However,  we  are 
unable  to  capture  the  exact  time  that  the  event  occurred.  Notice  that  the 
image  blurring  begins  at  t  =  0.1344sec  and  stops  only  after  t  =  0.201 6sec. 
This  is  caused  by  the  overlapping  of  the  edges  of  the  window  functions  at 
the  times  when  the  window  is  centered  in  the  period  that  the  blast  occurs. 
The  time  resolution  can  be  improved  with  shorter  window  width  or  sharper 
cutoff  at  the  edges  but  such  improvement  is  at  the  cost  of  poorer 
frequency  resolution. 

These  two  points  shall  be  discussed  in  more  detail  in  the  later  section. 

We  verify  the  consistency  of  the  result  obtained  with  a  second  case  in 
which  the  orientation  of  the  aircraft  is  set  to  face  away  from  the  radar  and  only  its 
back  area  is  subjected  to  the  blast  wave.  Consequently,  the  magnitude  of  the 
radial  velocity  variation  is  reduced  from  5.7m/s  to  less  than  0.7m/s.  (But  the  blast 
characteristic  remains  the  same.) 
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Figure  25  Acceleration  and  Velocity  variation  for  2nd  case 
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ISAR  image  from  FFT  for  2nd  Case 
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Figure  26 
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(b)  At  time  t  =  0.1536  to  0.3072sec. 


Figure  27  STFT  ISAR  images  for  2nd  case. 


Similar  cell  migration  effects  are  observed  in  this  case,  except  that  the 
extent  of  the  shift  is  slightly  lessened  compared  to  the  first  case.  This  behavior  is 
expected  since  the  velocity  variation  is  lowered  in  the  second  case.  The  cell 
range-migration  is  hardly  noticeable  (virtually  zero)  whereas  the  cross-range  cell 
variation  is  about  6-8  cells.  Also,  the  image  blurring  is  less  significant  in 
magnitude  and  duration  because  of  the  smaller  blast  force  and  shorter  blast 
duration.  Other  cases  were  also  simulated  and  examined  with  different  target 
orientations.  Their  results,  and  the  associated  conclusions,  do  not  differ  from 
those  described  above. 
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C.  SELECTION  OF  WINDOW  FUNCTIONS 

In  the  above  examples  we  have  seen  the  important  relationship  between 
the  time-domain  and  frequency-domain  resolutions  (a  relationship  which, 
unfortunately,  is  directly  linked  to  the  image  quality).  For  STFT,  these  resolutions 
are  determined  by  the  size  and  shape  of  the  window  function,  and  the  uncertainty 
principle  associated  with  the  Fourier  transform. 


1.  Resolution  of  Rectangular  and  Gaussian  Window  Functions 

In  general,  there  are  numerous  window  functions  to  choose  from  but,  for 
illustration  purposes,  we  will  restrict  our  attention  to  the  two  cases  for  linear 
TFT’s:  the  simple  rectangular  function;  and  the  Gaussian  function.  The 
rectangular  function  can  be  expressed  as 


w{t)  = 


,  a  a 
where  —  <  t  <  — 
2  2 


|1 

0  otherwise 

W (/)  =  a  sin  c{naf ) 


(6.1) 


where  a  is  the  width  of  the  window.  And  the  Gaussian  function  expression  is 
given  by 


Mt)  = 


cr2/2 

W(f)  =  (2a)I/2tr>/4e  2 


(6.2) 


The  resolutions  for  the  two  window  functions  —  in  both  time  and 
frequency  —  can  be  obtained  using  (Chen  [2]) 


At  = 


J  (t  ~  ju,)1  \H0f  dt 

—00 _ 

00 
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(6.3) 
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M==J - - 

I  Kof  *  jinoiv 

—00  -00 

Thus,  putting  (6.1)  into  the  Equations  in  (6.3),  we  derive  the  time  and 
frequency  resolutions  for  the  rectangular  function  to  be  At  =  a,  Af  =  1/a,  and  its 
time-bandwidth  product  AtAf  =  1.  Both  resolutions  are  directly  related  to  the 
width  of  the  window  and  multiplicative  inverses  of  each  other. 

Also,  putting  (6.2)  into  (6.3),  we  obtain  for  the  Gaussian  window  function, 
At  =  <j/^2,  Af  =  1  /  (<j^2)  and  AtAf  =  V2.  The  two  resolutions  are  related  to  the 
shape  of  the  window  function,  which  is  determined  by  a.  This  gives  us  an  added 
advantage  of  stretching  the  window  width  to  cover  the  entire  time-histories  of  M 
cells. 


2.  Effect  of  Using  Rectangular  Window  Function  in  STFT  Imaging 

In  our  earlier  demonstration  (the  first  case  shown  in  Figure  24),  we  used  a 
rectangular  window  function  with  a  window  size  a  =  ‘A  of  64  cells  =  16  cells.  The 
interval  between  each  cell  was  N  x  T2  =  64  x  150jusec  =  9.6msec.  By  setting  the 
window  width  to  16  cells,  we  obtained  the  resolutions  At  =  a  =  153.6msec  and  Af  = 
1/a  =  6.51  Hz. 

For  that  example,  the  blast  event  happened  between  time  steps  0.165  and 

0.19sec  and  was  only  25msec  in  duration.  Thus,  the  selected  time  resolution  is 

about  six  times  larger  than  the  blast  duration.  This  fact  explains  why  the  STFT 

ISAR  image  could  not  faithfully  capture  the  event.  The  reconstructed  images 

began  to  mark  the  blurring  effect  at  t  =  0.096sec  and  end  at  about  t  =  0.2496sec. 

The  difference  of  153.6msec  is  equivalent  to  the  time  resolution  of  the  rectangular 

function.  To  obtain  accurate  tracking  of  the  event,  the  window  size  should  be 

smaller  than  25msec,  which  (for  our  case)  is  one  to  two  m-cells.  But  in  doing  so, 

the  frequency  resolution  will  decline  significantly  to  about  1/ a  =  1/  (2  x  9.6msec)  = 
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52.1Hz.  At  this  frequency  resolution,  it  is  virtually  impossible  to  map  out  the  image 
at  all.  Below  is  an  example  showing  the  serious  degradation  of  the  image  by 
setting  the  window  size  to  2  m-cells. 
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Figure  28  STFT  ISAR  Time-history  images  with  width  set  to  2  m-cell. 


The  Doppler  resolution  for  image  quality  is  inversely  proportional  to  the 

integration  time  T.  Since  T  =  0.6144sec,  we  obtain  the  required  frequency 

resolution  to  be  AfD  =  1  /  T  =  1.62Hz.  But  the  frequency  resolution  of  the  above  is 

four  times  coarser  than  AfD.  As  a  result,  the  image  quality  is  generally  poorer 

compared  to  the  conventional  FFT  ISAR  imaging.  We  can  improve  the  image 

quality  by  extending  the  window  size  to  !4  of  the  M  cells  (or  more)  but  that  would 

mean  losing  the  time  resolution  and  the  dynamic  behavior  of  the  aircraft 

response  could  no  longer  be  tracked.  An  example  of  this  effect,  with  the  window 

size  set  to  55  m-cells,  is  shown  in  Figure  29.  Notice  that  the  position  offset  is 
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shown  right  at  the  start  of  the  integration  time  while  the  blast  effect  occur  at  t  = 
0.165sec.  But  the  quality  of  the  images  is  better  than  previous  ones. 
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Figure  29  STFT  ISAR  Time-history  images  with  width  set  to  55  m-cell. 


3.  Effect  of  Using  Gaussian  Window  Function  in  STFT  Imaging 

Compared  to  the  rectangular  window  function,  the  time-bandwidth  product 
of  the  Gaussian  window  of  only  >2  allows  better  time  and  frequency  resolution  to 
be  obtained  with  lower  limits  on  each.  Figure  30  shows  a  graph  comparing  the 
Gaussian  and  rectangular  window  function  in  term  of  their  time  versus  frequency 
resolution.  The  graph  is  obtained  using  their  respective  time-bandwidth  product. 
From  the  graph,  we  can  easily  see  that  the  Gaussian  window  requires  a  smaller 
window  size  to  obtain  quality  image  frequency  resolution.  Conversely,  for  the 
same  time  resolution  it  has  better  frequency  resolution  than  the  rectangular 
function. 
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Figure  30  Time  and  frequency  resolution  comparison  chart  for  rectangular  and 

Gaussian  window  functions. 
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(a)  At  time  t  =  0to  0.1536sec. 
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(b)  At  time  t  =  0.1536  to  0.3072sec. 


Figure  31  STFT  ISAR  Time-history  images  with  Gaussian  window  function. 


Figure  31  shows  the  same  32  STFT  images  generated  using  the 
Gaussian  window  and  should  be  contrasted  with  the  rectangular  window  used 
earlier.  To  enable  a  direct  comparison,  we  set  the  time  resolution  of  the  Gaussian 
window  to  be  the  same  as  in  the  first  (rectangular  window)  case  i.e.  16  x  9.6msec 
=  153.6msec.  Since  the  time  and  frequency  resolutions  are  determined  by  the 
shape  of  the  window  function  rather  than  its  width,  we  use  a=  At ^2  =  0.2172  to 
obtain  the  corresponding  frequency  resolution  as  Af  =  1  /  (a^/2)  =3.26Hz.  This 
value  is  only  half  of  the  previous  case. 
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We  can  see  that  the  time  resolution  remains  the  same  by  observing  the 
blurring  effect  starting  at  t  =  0.096sec  and  ended  slightly  early  at  t  =  0.2208sec. 
The  difference  of  about  124.8sec  is  close  to  the  time  resolution.  However,  the 
images’  quality  has  improved  over  that  of  the  earlier  case  (as  a  result  of  the  finer 
frequency  resolution). 


D.  EFFECT  OF  RADIAL  VELOCITY  VARIATION  ON  ISAR  IMAGE 

In  the  images  obtained  using  FFT  or  STFT,  cell  migration  in  both  range 
and  cross-range  directions  was  quite  noticeable  —  particularly  along  the  cross¬ 
range  axis.  This  migration  caused  the  ‘smearing’  of  the  image  and  the  blurring 
effect. 

To  quantify  the  extent  of  the  cell  migration  along  the  range  due  to  the 
presence  of  the  radial  velocity  component,  we  use  the  expression  (2.27)  given  in 
Chapter  II  to  obtain  a  plot  of  AN  versus  the  radial  velocity  (Figure  32).  For  a 
change  in  the  radial  velocity  that  corresponds  to  about  4-5  m/s  in  our 
demonstrated  blast  characteristic,  the  estimated  shift  from  Figure  32  below  is  AN 
=  4  range  cell.  This  is  about  the  same  amount  of  range  cell  shift  counted  in  the 
ISAR  images  (which  is  3  cells  given  in  the  previous  section).  Surprisingly,  it  is 
not  a  significant  contribution  to  image  distortion  but  only  causes  a  slight 
‘smearing’  along  the  range  cell  seen  earlier. 


Figure  32  Extent  of  cell  migration  due  to  velocity. 
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Comparatively,  cell  migration  along  the  cross-range  is  much  higher 
although  there  is  no  angular  variation  components  present  in  the  demonstration. 
To  explain  this  strange  effect,  consider  the  radar  signal  received  within  a 
particular  range  cell  index  n  given  below  (Chen  [2]). 

=  +  xcos0(t)- yksin0(t)]^  (6.4) 

where  Mn  is  the  total  number  of  scattering  points  that  share  the  same  range,  the 
target  aspect  angle  6(t)  is  expressed  similarly  to  (2.3),  and  R(t)  is  given  by  (5.8) 
as  shown  below. 

8(t)  =  0O  +  cot  +  ^yt2  + ....  (6.5) 


m=R0+v, 


sin  0B 
Mg 


\P(t',Rc)dt' 


(6.6) 


where  R0  is  the  range  at  the  start  of  data  collection  (integration  period)  and  Vt  is 
the  initial  radial  velocity.  From  Chapter  II,  the  cross-range  cell  migration  is 
maximal  when  the  initial  angle  is  zero.  Setting  60  =  0  and  assuming  that  there  is 
no  angular  acceleration  yields  6(t)  =  cot.  Substitute  (6.6)  into  (6.4)  and,  for  zero 
initial  velocity  (as  in  our  demonstrated  case),  we  get 
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By  expanding  sinntf  =  ®r-£yV/ 6  +  .. .and  coscot  =  \-ort2 1 2  + ...  in  (6.7)  and 
ignoring  the  higher  order  terms,  we  obtain 
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It  is  easier  to  quantify  the  effect  on  cross-range  cell  migration  if  we 
consider  the  case  where  there  is  only  one  scattering  point.  The  Doppler  shift  in 
the  range  cell  with  index  n  then  can  be  obtained  by  differentiating  the  range 
related  phase  terms  in  (6.8)  to  give 


where 


fa,  =  4A(0 

at 


dt  c 


An  sin  0R  c 
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(6.9) 


for  t<ta 

for  ta  <t  <td  +ta  (6.10) 
for  td+ta<t 


and  p0,  td,  ta  and  a  depend  on  Rc,  the  distance  between  the  blast  point  and  the 
center  of  mass  of  the  aircraft.  Note  that  the  second  term  Ap  sin  0B  P(t,  Rc)  /  MG  in 
equation  (6.9)  represents  the  impulse  response  of  the  aircraft,  and  since  the  blast 
duration  is  much  shorter  than  the  integration  period,  this  term  is  independent  of 
T.  Using  Equation  (2.40)  with  Arc  =  c  /  2fccoT,  we  can  estimate  the  amount  of 
cross-range  cell  shift  to  be 
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where  I/A  can  be  obtained  using  (3.8)  for  known  pQ,  td  and  a.  From  the  simulation 
parameters  given  in  the  program  in  Appendix  A  —  pa  =  556274. 3N/m2 ,  td  = 

0.01 41  sec,  a  =  2.53,  we  obtain  I/A  =  1974.7N-sec/m2 . 

Also,  in  the  first  demonstration  case,  we  have  Arc  =  1.23m,  Ap  =  64m2,  0B  = 
-90,  co  =  3. 7 7 deg/ sec  =  0. 06 5 8r ad/sec,  MG  =  22000kg.  Using  Equation  (6.1 1),  we  get 
AM  =  72  cells.  This  is  more  than  the  total  number  of  cells  in  the  cross-range. 
Hence,  the  scattering  points  will  crossover  into  the  subsequent  time-history 
series  of  the  image  and  ‘appear’  as  if  they  belong  the  scattering  points  in  that 
time-history.  In  this  case,  AM  -  M  =  72  -  64  =  8  cells,  which  is  about  the  same 
number  of  cells  that  was  observed  in  the  first  demonstrated  case.  From  Equation 
(6.11),  the  extent  of  the  cross-range  cell  migration  seems  to  depend  on  the 
aircraft  impulse  response  to  the  blast  wave.  This  has  been  verified  also  for  other 
blast  orientation  and  lower  blast  magnitude. 
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VII.  CONCLUSION  AND  RECOMMENDATIONS 


A.  CONCLUSION 

The  work  presented  in  this  thesis  is  a  preliminary  demonstration  of  the 
possibility  of  using  ISAR  to  conduct  battle  damage  assessment  (BDA)  (based  on 
computer  simulation).  The  significant  advantage  of  ISAR  —  as  opposed  to  other 
types  of  sensors  —  is  its  ability  to  produce  high  resolution  target  images  at 
greater  surveillance  ranges  than  optical  systems.  With  high  quality  images,  the 
positive  evaluation  of  the  extent  of  the  target  damage  can  be  quickly  made  and 
immediate  follow-on  action  can  be  taken  if  required. 

However,  most  ISAR  images  are  generally  static  in  nature  and  require 
relatively  long  processing  periods  to  construct.  The  quality  of  the  image  is  often 
tied  to  the  processing  time  and  the  prior  conditions  associated  with  the 
processing.  The  foremost  of  these  conditions  is  the  assumption  that  the  received 
radar  signal  is  reflected  from  a  target  with  only  constant  angular  velocity 
components  (an  assumption  that  is  used  to  map  the  Doppler  variation  of  the 
respective  scattering  points  caused  by  the  angular  rotation  to  the  cross-range 
position  of  the  scattering  points).  The  presence  of  the  radial  velocities  and 
deviation  of  the  angular  rotation  rate  will  lead  to  image  distortion.  Unfortunately, 
these  variations  cannot  be  ignored  during  BDA. 

In  this  study,  our  main  interest  is  on  aircraft.  When  an  aircraft  is  hit  by 
explosive  munitions  two  distinct  behaviors  can  be  observed.  First,  the  blast  wave 
of  the  explosion  will  interact  with  the  surface  body  of  the  aircraft  and  cause  the 
aircraft  to  have  translational  motion,  rotational  motion  and  vibration  responses. 
The  characteristics  of  these  responses  depend  on  the  magnitude  of  the  blast 
loading,  the  orientation  of  the  aircraft  surface  to  the  blast  wave  and  the  direction 
of  these  responses  with  respect  to  the  ISAR  observing  the  target.  Second,  most 
explosive  munitions  against  airborne  targets  carry  pre-fabricated  fragments  and 
shrapnel  aimed  to  penetrate  the  thin  layer  of  the  target  to  cause  a  variety  of 
damages:  from  loss  of  flight  control,  tearing  of  structures  to  fuel  tanks  and 
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munitions  detonation.  The  nature  of  this  latter  effect  is  more  difficult  to  account 
for  in  simulation. 

1.  Approach  to  Model  the  Aircraft  Response  When  Hit 

Our  approach  was  to  examine  the  response  of  the  aircraft  in  a 
deterministic  manner  by  looking  at  blast  loading  effect.  Even  so,  there  was  a 
need  to  narrow  down  the  scope  of  the  study  to  within  the  time  frame  of  this  thesis 
work  and  our  main  focus  was  solely  on  the  translational  and  rotational  motions  of 
the  aircraft.  This  is  a  reasonable  approach  if  we  consider  the  aircraft  as  a  rigid 
body  in  air.  These  motions  are  likely  to  be  the  dominant  responses  apart  for  the 
latter  effect  caused  by  the  fragment  penetrations.  Besides,  they  are  directly 
detectable  by  the  ISAR  system. 

Our  simulation  begins  with  the  formation  of  a  blast  wave  in  air  and 
requires  its  characteristic  parameters  to  be  understood.  They  are  discussed  in 
the  initial  sections  of  Chapter  III.  It  was  shown  that  the  blast  wave  can  be 
represented  by  an  impulse  pressure  force  as  a  function  of  time  by  Equation  (3.1) 
that  depends  on  four  key  parameters  —  its  peak  overpressure,  blast  duration, 
arrival  time  and  its  waveform  parameter.  These  four  parameters  are,  in  turn, 
dependent  on  the  distance  between  the  center  of  blast  and  the  point  of 
interaction  between  the  pressure  force  and  the  surface.  The  parameters  are 
usually  described  using  empirical  equations. 

The  overall  translational  and  rotational  responses  of  the  aircraft  from  the 
blast  loading  can  be  derived  using  the  equation  of  motions  (Newton’s  2nd  law). 
However,  the  non-symmetrical  geometry  of  the  aircraft  surface  and  the  unequal 
distribution  of  the  blast  wave  (as  a  result  of  its  spherical  divergence)  increase  the 
complexity  of  the  problem  and  it  has  been  shown  in  Chapter  III  that  it  is  difficult  to 
yield  an  analytical  closed-form  solution.  However,  the  results  of  Equations  (3.21) 
and  (3.22)  can  be  obtained  through  the  use  of  extensive  numerical  simulation. 
The  methodological  approach  for  the  simulation  was  also  briefly  outlined  in  the 
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chapter.  Due  to  the  time  constraint  of  this  thesis  work,  we  were  (regrettably)  not 
able  to  realize  it. 

Instead,  we  simplified  the  study  by  modeling  the  translational  response  of 
the  aircraft  alone.  This  restriction  can  be  accomplished  with  reasonable 
assumptions  made  concerning  the  nature  of  the  blast  wave  and  the  surface  area. 
We  have  assumed  a  planar  blast  wave  with  equal  pressure  distribution  across 
the  aircraft  body.  Also,  the  surface  area  of  the  aircraft  during  the  evaluation 
period  does  is  assumed  to  not  change  with  time.  These  assumptions  help  to 
decouple  the  time  and  geometrical  surface  dependent  pressure  distribution  of 
(3.21)  into  two  separable  variables  shown  in  Chapter  IV  and  in  (5.1)  —  both 
independent  of  each  other.  The  translational  motion  can  be  obtained  in  terms  of 
radial  acceleration,  velocity  and  its  resultant  range  variation  with  respect  to  time, 
by  solving  (5.1)  analytically.  However,  it  is  not  possible  to  similarly  obtain  the 
angular  velocity  variation  since  the  resultant  changes  in  angular  momentum 
depend  on  the  differences  in  pressure  distribution  across  the  surface. 


2.  Approach  to  Overcome  the  Image  Distortion 

The  blast  loading  on  the  aircraft  induces  variations  in  both  the  radial  and 
angular  velocities.  The  velocity  variations  are  encoded  as  changes  in  the 
Doppler  frequencies  of  each  returned  scattering  signal.  This  information  is 
captured  by  the  ISAR  and  use  for  radar  imaging.  However,  the  time  resolution 
(which  is  equivalent  to  its  integration  time)  of  the  conventional  Fourier  method  of 
image  processing  is  extremely  poor  and  it  is  difficult  to  extract  the  dynamic 
frequency  changes  during  and  after  the  blast  event  occurs.  As  a  result,  the 
image  becomes  distorted  due  to  the  inability  of  the  Fourier  processing  to  account 
for  the  velocity  variations. 

On  the  other  hand,  the  Time-frequency  Transform  (TFT)  method  of  signal 
analysis  offers  an  added  dimension  allowing  the  signal  to  be  represented  as  its 
frequency  changes  with  time.  We  have  shown  through  the  use  of  a  simple  TFT 
method  of  analysis  (which  is  the  Short-time  Frequency  Transform  or  STFT)  that 
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the  dynamic  behavior  of  the  aircraft  as  it  is  being  hit  can  be  revealed  in  the  ISAR 
images.  The  TFT  images  provide  more  information  about  the  event  than  with 
FFT  processing  alone.  The  use  of  TFT  for  ISAR  imaging  is  also  beneficial  to 
similar  future  work  because  of  its  ability  to  extract  time-varying  frequencies  since 
many  of  these  efforts  will,  in  some  way  or  another,  involve  examining  Doppler 
variations  of  measured  radar  signals. 


B.  RECOMMENDATIONS  FOR  FURTHER  WORK 

Although,  the  use  of  ISAR  systems  with  TFT  processing  has  been 
demonstrated  in  this  thesis  for  BDA,  the  simulation  model  used  for  proof-of- 
concept  gives  only  an  initial  indication  of  this  possibility.  Further  work  is  highly 
recommended  to  refine  the  present  study  in  the  three  major  areas  discussed 
below. 


1.  Complete  Model  of  the  Blast  Loading  on  the  Aircraft 

In  this  simulation,  only  the  simplest  case  —  the  translational  motion  of  the 
aircraft  —  is  accounted  for.  Even  then,  there  were  certain  assumptions  made  to 
reduce  the  magnitude  of  the  problem  so  that  some  rough  estimates  could  be 
obtained  (for  example,  the  possible  angular  variations  were  not  considered).  To 
enable  a  complete  description  of  the  aircraft  response  to  the  blast,  there  is  a 
need  to  develop  a  full  model  that  accounts  for  all  the  parameters  given  in 
Equations  (3.21)  and  (3.22).  It  is  not  an  impossible  task,  but  requires  careful 
consideration  of  the  geometry  of  the  problem  between  the  aircraft  (its  surface 
distribution  and  orientation),  the  location  of  the  blast  point  and  its  pressure 
distribution  in  time  and  space,  and  the  position  of  the  observing  radar.  The 
problem  can  also  be  extended  from  the  present  two-dimensional  analysis  to  a 
real  world  three-dimensional  case  (but  also  with  an  increasing  level  of 
complexity).  As  discussed  in  Chapter  III,  this  blast  effect  model  can  be  developed 
as  a  separate  model  and  its  results  can  be  simply  expressed  as  a  collection  of 
time  varying  range  and  angular  information  over  the  entire  integration  period 

(following  the  methodological  approach  proposed). 
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2.  Study  of  Target  Break  Up  Process 

A  further  extension  of  the  blast  model  is  to  include  the  target  break  up 
process  due  to  the  penetration  of  the  fragments  into  the  aircraft.  However,  the 
break  up  process  is  generally  treated  statistically  and  depends  on  the  details  of 
the  fragment  penetration.  As  a  crude  start,  we  could  divide  the  aircraft  into 
different  sections  and  each  section  could  be  driven  by  different  translational  and 
rotational  responses  induced  by  the  blast  effect  due  to  the  difference  in  pressure 
wave  distribution  in  each  section.  In  this  case,  a  series  of  range  and  angular 
variations  for  different  sections  can  serve  as  input  to  the  ISAR  imaging  model. 
The  ISAR  imaging  model  will  then  introduce  the  appropriate  phase  variations 
according  to  the  range  and  angular  inputs  to  respective  scatterers  categorized  by 
the  different  sections. 

3.  Improving  the  Time  and  Frequency  Resolution  of  TFT 

While  the  linear  Short-time  Frequency  Transform  (STFT)  is  the  simplest  to 
implement  for  our  demonstration,  we  have  seen  that  the  time  and  frequency 
resolutions  still  lack  the  sensitivity  required  to  provide  faithful  tracking  of  the 
aircraft  dynamic  behavior  due  to  blast  and  simultaneously  provide  high  quality 
radar  imagery.  Depending  on  the  choice  of  the  sliding  window  width  and  shape, 
the  time-resolution  can  be  adjusted  in  such  a  way  that  the  dynamic  motion  of  the 
target  within  the  integration  period  can  be  extracted  and  tracked.  But  a  major 
drawback  of  this  approach  is  that,  by  gaining  time  resolution,  the  frequency 
resolution  will  necessarily  suffer  —  and  the  quality  of  the  image  depends  on  the 
frequency  resolution.  This  limitation  is  fundamental  and  is  imposed  by  the 
uncertainty  principle  of  time  and  frequency  in  any  Fourier  process.  For  a  linear 
TFT,  the  best  resolution  is  achievable  using  Gaussian  window  with  the  time- 
bandwidth  product  of  V2. 

From  a  system  point  of  view,  one  can  further  improve  both  resolutions  by 
having  a  greater  number  of  range  and  cross-range  cells  to  represent  the  image, 

although  other  radar  considerations  (such  as  the  carrier  frequency,  pulse  interval 
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and  frequency  step  size  with  respect  to  the  size  of  the  target  and  its  rotation  rate) 
will  have  to  be  factored  in  as  they  are  discussed  in  Chapter  II. 

The  alternate  approach  is  to  work  on  enhancing  TFT  image  processing. 
Most  efforts  involving  time-frequency  analysis  utilize  the  more  complex  variations 
of  STFT  to  overcome  the  time-frequency  resolution  limit  imposed  by  the  STFT. 
Examples  of  some  of  these  TFT  schemes  include  the  Continuous  Wavelet 
Transform  (CWT),  the  bilinear  Wigner-Ville  Distribution  (WVD),  the  Cohen’s  class 
series  and  adaptive  time-frequency  representations.  Many  of  these  TFT 
approaches  have  reported  better  time-bandwidth  products  and,  of  these,  the 
WVD  seems  to  give  best  time-frequency  resolution  (but  at  the  expense  of 
significant  cross  interference  terms). 
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APPENDIX  A:  MATLAB  CODES 


%  TFT  ISAR  IMAGE  SIMULATOR 

%  This  program  is  developed  to  demonstrate  the  possibility  of  reconstructing 
%  ISAR  image  using  Time-Frequency  Transform  Method.  It  can  be  used  for 
%  viewing  the  dynamic  behaviour  of  the  target  due  to  blast  effect. 

%  It  contains  4  parts: 

%  a.  Part  1  -  Define  the  radar  parameters  and  Target  Model. 

%  b.  Part  2  -  Generate  the  blast  effect  on  traslational  motion. 

%  c.  Part  3  -  Creates  the  returned  stepped  frequency  signal. 

%  (This  is  modified  from  the  original  code  contained  in  Jae  S.  S.  et 

%  al  "Range-Doppler  Radar  Imaging  and  Motion  Compensation",  Norwood,  MA, 

%  Artech  House,  2001). 

%  b.  Part  4  -  Construct  the  ISAR  Image  Using  2D  FFT  and  TFT. 

%  Assumptions  Made:  Range  aligned.  Motion  compensated.  Zero  receiving  system  delay 
%  and  constant  amplitude  scatterers  (unity). 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Clear  Memory 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
clear  all;  fig_con  =  0; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  PART  1:  Define  the  radar  parameters  and  Target  Model. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Define  radar  parameters 


M  =  64; 

N  =  64; 
stepf  =  2.2e6; 
fc  =  3e9; 

Ro  =  22567.2; 
Vro  =  0; 
Accra  =  0; 
Tetao  =  0; 
wo  =  3.77; 
alo  =  -0.03; 
c  =  3e8; 

T2  =  150e-6; 
T1  =  l/step_f; 


%  Total  number  of  bursts 
%  Total  number  of  pulses 
%  Frequency  step  (Hz) 

%  Starting  frequency  (Hz) 

%  Initial  range  (meters) 

%  Initial  translational  velocity  (m/ses) 

%  Initial  translational  aceleration  (m/sec-sq) 
%  Initial  angle  of  rotation  (deg) 

%  Initial  angular  velocity  (deg/sec) 

%  Initial  angular  acceleration  (deg/sec-sq) 

%  Speed  of  EM  wave  (m/sec) 

%  Pulse  repetition  interval  (sec) 

%  Pulse  length  (sec) 


wo  =  wo*pi/180;  alo  =  alo*pi/180;  Teta_o  =  Teta_o*pi/180; 

T1  =  l/step_f; 

%  Define  two  dimensional  density  reflectivity  function  of  Aircraft 
%  Rows  are  x  and  y  coordinates, 

TgtModel  =  zeros(N,M); 

Tgt_Model(29, 14)=  1 Tgt_Model(3 1 , 1 4)=1 ; 
Tgt_Model(27,16)=l;,Tgt_Model(33,16)=l; 

Tgt_Model(27, 1 9)=1 Tgt_Model(3 3 , 1 9)=1 ; 
Tgt_Model(27,23)=l;,Tgt_Model(33,23)=l; 
Tgt_Model(25,25)=l;,Tgt_Model(35,25)=l; 
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Tgt_Model(23,27)=l; 

Tgt_Model(21,29)=l; 

Tgt_Model(19,31)=l; 

Tgt_Model(17,33)=l; 

Tgt_Model(15,35)=l; 

Tgt_Model(13,37)=l; 

Tgt_Model(15,38)=l; 

Tgt_Model(17,39)=l; 

Tgt_Model(19,38)=l; 

Tgt_Model(21,37)=l; 

Tgt_Model(23,36)=l; 

Tgt_Model(25,35)=l; 

Tgt_Model(27,34)=l; 

Tgt_Model(27,37)=l; 

Tgt_Model(27,40)=l; 

Tgt_Model(28,43)=l; 

Tgt_Model(26,45)=l; 

Tgt_Model(24,47)=l; 

Tgt_Model(28,47)=l; 

Tgt_Model(30,42)=l; 


Tgt_Model(37,27)=l; 

Tgt_Model(39,29)=l; 

Tgt_Model(41,31)=l; 

Tgt_Model(43,33)=l; 

Tgt_Model(45,35)=l; 

Tgt_Model(47,37)=l; 

Tgt_Model(45,38)=l; 

Tgt_Model(43,39)=l; 

Tgt_Model(41,38)=l; 

Tgt_Model(39,37)=l; 

Tgt_Model(37,36)=l; 

Tgt_Model(35,35)=l; 

Tgt_Model(33,34)=l; 

Tgt_Model(33,37)=l; 

Tgt_Model(33,40)=l; 

Tgt_Model(32,43)=l; 

Tgt_Model(34,45)=l; 

Tgt_Model(36,47)=l; 

Tgt_Model(32,47)=l; 


figcon  =  figure(fig_con+l); 
imagesc(-Tgt_Model),  colormap(gray) 
ylabel('Range  Bin'),  xlabel('Cross-Range  Bin') 

title('Density  Reflectivity  Distribution  of  the  Target  Scattering  Points') 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Part  2:  Generate  the  blast  effect  on  traslational  motion. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  Aircraft  Area  Profde  -  Front  and  Side 
Areafront  =  7.6416;  %  Front  Area  in  m-sq 
Areaside  =  64;  %  Side  Area  in  m-sq 
Mass_Aircraft  =  22000;  %  Mass  of  Aircraft  in  kg 

%  Blast  Parameters  -  Air  Explosive  Profile. 

%  Based  on  Medium  Range  SAM  -  Explosive  Weight  =  35kg,  Miss  Distance  =  5m, 

%  Target  Fleight  =  3000m,  TNT  Equivalent  =1.5  (FIE) 
fract  =  0.25; 

Ar  Delay  =  fract*N*M*T2;  %  Arbitaiy  Delay  in  fraction  of  integration  time 
%  Obtain  from  Excel  Table  -  "Blast  Effect  Study  -  Input  Parameters.xls" 

Po  =  556274.3;  %  Peak  Overpressure  in  Pa  (N/m-sq) 

ta  =  0.015593  +  Ar_Delay;  %arrival  time  in  sec 
td  =  0.0141 16;  %blast  duration  in  sec 
alpha  =  2.53;  %  waveform  parameter 

%  Geometry  of  Blast  -  Equation  (5.10),  (5.7),  (5.8)  &  (5.9) 

TetaBdeg  =  -90;  %  direction  of  blast  (deg)  CW  ref  to  neg.  crossrange  axis 

TetaB  =  TetaB  deg  /  180  *  pi; 

AreaP  =  abs( Areafront  *  cos(TetaB  +  Teta  o))  +  abs(Areaside  *  sin(TetaB  +  Teta  o));  %  Projected  Area 
Fyo  =  Po  *  AreaP  *  sin(TetaB); 

Ayo  =  -  Fyo  /  Mass_Aircraft; 

%  Sampling  times  Si 
%  Sik  vector  of  64x64  elements 
i  =  1:M*N; 
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Sik  =  (Tl/2)  +  2*Ro/c  +  ((i-l)*T2);  %  Using  Equation  (2.9) 


%  Range  Variation  (  T  due  to  blast  effect,  'O'  linear  expression) 

WithBlast  =1;  %  Blast  Switch 

if  WithBlast  %  Using  Equation  (5.7),  (5.8)  and  (5.9) 

Rik  =  zeros(l,N*M); 
for  vect_index  =  1  :M*N 
if  (Sik(vect_index)  <=  ta) 

Rik(vectindex)  =  Ro  +  Vro*Sik(vect_index); 
elseif  (Sik(vect_index)  >  ta)  &&  (Sik(vect_index)  <=  (ta  +  td)) 

Rik(vect_index)  =  Ro  +  Vro*Sik(vect_index)  +  Ayo*((td/alpha)A2  *  (2/alpha-l)*(l-exp( 
-alpha*(Sik(vect_index)-ta)/td))  +  td/alphaA2  *  (Sik(vect_index)-ta)*( 
alpha-1-  exp(-alpha*(Sik(vect_index)-ta)/  td))); 

else 

Rik(vect_index)  =  Ro  +  Vro*Sik(vect_index)  +  Ayo*((td/alpha  -  td/alphaA2  *  (1  -  exp( 
-alpha)))*(Sik(vect_index)  -  (ta  +  td))  +  tdA2/alphaA3  *  (2  -  alpha)  *  ( 

1  -  exp(-alpha))  +  tdA2/alphaA2  *  (alpha  -  1  -  exp(-alpha))); 
end 
end 
else 

Rik  =  Ro  +  Vro*Sik  +  0.5*Acc_ro*Sik.A2;  %  Equation  (2.2) 
end 

Rik  =  reshape(Rik,N,M); 

%  Angular  positions  -  Equation  (2.3) 

Tetik  =  Teta_o  +  wo*Sik  +  0.5*alo*Sik.A2; 

Tetik  =  reshape(Tetik,N,M); 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Part  3 :  Creates  the  returned  stepped  frequency  signal. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Frequency  steps  -  Equation  (2.8) 
i  =  1:M; 

fil  =  fc  +  (i-l)*step_f;  fi  =  []; 
for  i  =  1:N  fi  =  [fi  fil];  end 
clear  fil; 

fis  =  reshape(fi,N,M); 

%  Create  the  target  image  by  summing  all  the  returns  together. 

%  Step  1:  Find  the  number  of  point  scatterers  on  the  image 
numpo  =  length(find(Tgt_Model  ==  1)); 

[xx  yy]  =  find(Tgt_Model  ==  1); 
pos  =  1; 

target  =  zeros(M,N); 

%  Step  2:  Perform  the  summation  -  Equation  (2.10) 
for  i  =  1  :numpo 

target  =  target  +  exp(-j*4*pi/c*fis.*Rik).*exp(-j*2*pi*(  (xx(pos)-N/2)*2*(fis/c).*cos(Tetik) 

-  (yy(pos)-N/2)*2*(fis/c).*sin(Tetik))); 
pos  =  pos  +1; 
end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  PART  4:  Construct  the  ISAR  Image  Using  2D  FFT  and  TFT. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  TFT  Cube  Initialisation 

stft  cube  =  zeros([N  M  N]);  %  Instant  Frequency-Time  Flistory-Range  Cell 
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%  Obtain  the  Synthetic  Range  Profile  of  the  target  -  Equation  (2.11) 
th_target  =  fft(target);  %  fft  on  column  of  matrix 

%  Construct  the  Conventional  FFT  Image  -  Equation  (2.12) 

figcon  =  figure(fig_con+l); 

imagesc(fftshift(abs(fft(th_target. ').'))); 

ylabel('Range  Bin') 

xlabel('Cross-Range  Bin') 

colormap(gray); 

%  Perform  TFT  on  the  doppler  history  of  each  Range  cell  index  using  STFT 
%  Ref:  Chapter  IV  Section  D 
%  Window  function 

h  =  window(45, 'gauss', 0.005);  %  For  gauss  -  need  to  adjust  width  according  to  sigma  /  0.005 
for  vectindex  =  1  :M 

stft_cube(:,:,vect_index)  =  tfrstft(th_target(vect_index,:).',l:M,M,h).'; 
end 

%  Generate  and  Plot  the  Reconstructed  ISAR 
for  vect_indexl  =  1:4 

fig_con  =  figure(fig_con+l); 
for  vect_index  =1:16 

fii  =  fftshift(squeeze(stft_cube(vect_index+(vect_indexl-l)*16,:,:)).'); 
subplot(4,4,vect_index),  imagesc(abs(fii)),  colormap(gray); 

TFlmarker  =  vect_index+(vect_indexl-l)*16; 

title(['TFl',num2str(TF[marker),'  (t  =  ',  num2str(TFlmarker*T2*N),'sec)']); 
axis  off; 
end 
end 
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